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METHODS AND APPARATUS FOR ANALYSING ULTRASOUND IMAGES 

FIELD AND BACKGROUND OF THE INVENTION 

The present invention relates to image analysis and, more particularly, to 

5 methods and apparatus for analyzing and improving the quality of images obtained, 
e.g., during echocardiography. 

Imaging is recognized as the most commonly used tool in medical diagnostics. 
Due to its non-invasive nature, imaging is the preferred procedure prior to any 
invasive treatment or analysis. Medical imaging techniques include Magnetic 

10 Resonance Imaging, X-ray imaging, gamma imaging, ultrasound imaging and the like. 
In Magnetic Resonance Imaging, magnetic fields interact with the spins of the atoms 
in a tissue and the interaction results are monitored and analyzed to provide an image 
of the tissue. In X-ray imaging, X-ray radiation is applied to the body and different 
absorption and transmission characteristic of different tissues generate an image 

15 thereof. In gamma imaging, a radioactive isotope is injected to, inhaled by or ingested 
by a patient. The isotope is chosen based on bio-kinetic properties that cause 
preferential uptake by different tissues. Radiation emitted by the radioactive isotope is 
detected by radiation detectors outside the body, giving its spatial uptake distribution 
within the body. In ultrasound imaging, high frequency pulsed and continuous sound 

20 waves are applied to the body and the reflected sound waves are used to develop 
images of internal organs and the vascular system. The sound waves are generated 
and recorded by transducers or probes that are either passed over or inserted into the 
body. The resulting images can be viewed immediately on a video display or can be 
recorded for later evaluation in a single image or a cine-loop format. 

25 Diagnostic ultrasound imaging is presently a preferred imaging modality in 

many medical fields, such as radiology, cardiology and gynecology. Cardiologists and 
other medical practitioners use cardiac ultrasound imaging, also termed 
echocardiography, to evaluate the condition of the heart. Echocardiography has the 
advantage of being a non-invasive, quick, inexpensive, convenient and safe diagnostic 

30 procedure, and is therefore practiced in many hospitals as well as private clinics. 

The primary drawback of echocardiography is the difficulty of acquiring good 
quality images in patients with poor acoustic windows. Moreover, clutter and poor 
resolution can compromise the clinical utility of images of any patient produced by 
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even the most sophisticated ultrasound scanners. With echocardiography, the 
difficulty of acquiring acceptable images is further compounded by the fact that the 
region of interest, the heart, has complex motion patterns. The advantages of 
echocardiography procedure on the one hand, and the unsatisfactory image quality on 

5 the other hand, have led researches to apply image processing techniques so as to at 
least partially improve the echocardiograph image. 

One image processing technique is "thresholding," in which one or more 
parameters (typically intensity thresholds), are used to generate an output image. For 
example, in one implementation of the thresholding procedure, the intensity of each 

10 pixel in the original image is compared with a single intensity threshold. The original 
image is mapped onto a binary image in which each pixel has one of two polarities 
(say "0" and "1") depending whether the intensity of the corresponding pixel in the 
original image is higher or lower than the intensity threshold. 

Image processing oftentimes involves mathematical operations performed on 

15 histograms characterizing the images. In image processing context, a histogram 
typically refers to a graph showing the number of pixels in an image at each different 
intensity value found in that image. For example, for an 8-bit grey-scale image, there 
are 2 8 = 256 different possible intensity values, and the histogram graphically displays 
256 numbers showing the distribution of pixels amongst those intensity values. 

20 A well-known mathematical operation is "histogram equalization" (see, e.g., 

Eltoft et al. "Real-Time Image Enhancement in Two-Dimensional Echocardiography", 
Computers in Cardiology 1984:481). This technique is based on the assumption that 
images embodying the maximal possible intensity range display optimal contrasts. In 
conventional histogram equalization, an intensity transformation, also known as a 

25 Brightness Transfer Function (BTF), is used to increase the spread of the intensity 
histogram characterizing the image. Histogram equalization can also be combined 
with thresholding and/or intensity transformation [Gonzalez R.C. and Woods R.W., 
"Digital Image Processing" Addison- Wesley, pp. 166-171, 1992]. However, the 
results obtained using the above techniques are far from being satisfactory. In 

30 particular, echocardiograph images processed using prior art techniques suffer from 
poor resolution and a substantial amount of noise. 

The major artifact in ultrasound images is clutter, which includes irrelevant 
information that appears in the imaging plane, obstructing the data of interest. There 
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are several causes for the appearance of clutter in an ultrasonic image. A first cause is 
effective imaging of off-axis objects, primarily due to highly reflective objects in the 
transducer's sidelobes (e.g., the ribcage and the lungs). A second cause is known as 
multi-path or reverberations. Due to the geometry of the scanned tissue with respect 
5 to the transducer, and the local reflective characteristics of the tissue, a substantial 
amount of the transmitted energy is bounced back and forth in the tissue before 
reaching the transducer. As a result, the signal measured at a specific range-gate 
includes contributions from incorrect ranges, in addition to the relevant range swath. 

A known method for handling clutter, in particular in ultrasound images of 

10 patients with low echogenicity is "contrast echocardiography" [Krishna et al, 
"Subharmonic Generation from Ultrasonic Contrast Agents," Physics in Medicine and 
Biology, 44:681, 1999]. In contrast echocardiography the backscatter from blood is 
enhanced to improve its delineation from the surrounding tissue. Another imaging 
method known to reduce clutter is harmonic imaging [Spencer et al, "Use of 

1 5 Harmonic Imaging without Echocardiographic Contrast to Improve Two -Dimensional 
Image Quality," American Journal of Cardiology, 82:794, 1998]. In harmonic 
imaging, the ultrasound waves are transmitted at one frequency and receiving at twice 
the transmitted frequency. These techniques however provide less than optimal 
results. Additionally, being based on adapting the data acquisition process, these 

20 techniques cannot be applied to all types of echocardiograph images. 

Several clutter rejection algorithms have been specifically developed for color- 
Doppler flow images in which effects of slow-moving objects are suppressed 
assuming that the blood flow velocity is much higher than the motion velocity of the 
surrounding tissue (to this end see, e.g., Herment et al, "Improved Estimation of Low 

25 Velocities in Color Doppler Imaging by Adapting the Mean Frequency Estimator to 
the Clutter Rejection Filter," IEEE Transactions on Biomedical Engineering, 43:919, 
1996; Bjaerum et al, "Clutter filters adapted to tissue motion in ultrasound color flow 
imaging," IEEE Transactions on Ultrasonics Ferroelectrics & Frequency Control, 49, 
6:693, 2002; Cloutier et al, "A new clutter rejection algorithm for Doppler 

30 ultrasound," IEEE Transactions on Medical Imaging, 22, 4:530, 2003; and Yoo et al, 
"Adaptive Clutter Filtering for Ultrasound Color Flow Imaging," Ultrasound in 
Medicine and Biology, 29, 9:131 1, 2003). 



It is recognized that the effectiveness in diagnostic imaging depends on the 
ability to accurately recognize the imaged organs. For example, in echocardiography, 
the determination of the location of the cardiac muscle within the scanned plane, and 
specifically of the Left Ventricle (LV), is of great importance. Information about the 

5 LV outlines as a function of time enables automatic extraction of rich local 
quantitative functional information. 

However, with the present signal-to-noise ratio and substantial amount of 
clutter in echocardio graph images, visual as well as automatic determination of the LV 
outlines is rather difficult. An inherent problem with automatic determination of the 

10 LV outlines is the complex motion of the Mitral Valve and the Papillary Muscles, 
which further increase the computational load. An additional problem is the 
significant variations between different patients and different measurements of the 
same patient. Several attempts have been made to develop algorithms for automatic 
detection of the LV outlines [US Patent Nos. 5,457,754 and 6,346,124; and Jacobs et 

15 al. , "Evaluating a Robust Contour Tracker on Echocardiographic Sequences," Medical 
Image Analysis, 3:63, 1999]. These attempts, however, had only limited success in 
border detection. For example, prior art fail to accurately detected contours outlining 
the outer boundaries of the LV. 

Other prior art of interest include, Ohyama et at., "Automatic Left Ventricular 

20 Endocardium Detection in Echocardiograms Based on Ternary Thresholding Method", 
IEEE Proceedings of 15th International Conference on Pattern Recognition 4:320-323, 
2000; and Abiko Y and Ito T, Nakajima M., "Improvement on Quality of 
Echocardiograms", Acoustical Imaging 23:169-176, 1997. 

It will be appreciated that there is a widely recognized need for, and it would 

25 be highly advantageous to have methods and apparatus for analyzing and improving 
the quality of images, devoid of the above limitations. 

SUMMARY OF THE INVENTION 

According to one aspect of the present invention there is provided a method of 
30 improving an image by transforming an intensity histogram thereof, the method 
comprising: (a) fitting the intensity histogram to a sum of a plurality of localized 
functions; (b) using the plurality of localized functions to define a plurality of 
localized intensity histograms; (c) for each localized intensity histogram, performing 
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at least one image enhancement procedure, thereby providing a plurality of improved 
localized intensity histograms; and (d) combining the plurality of improved localized 
intensity histograms, thereby transforming the intensity histogram of the image. 

According to further features in preferred embodiments of the invention 
described below, the method further comprises removing clutter from the image. 

According to still further features in the described preferred embodiments the 
plurality of localized functions comprises a first localized function, a second localized 
function and a third localized function, and further wherein the plurality of localized 
histograms comprises a first localized histogram, a second localized histogram and a 
third localized histogram. 

According to still further features in the described preferred embodiments the 
method further comprises reducing a range of the first localized histogram. 

According to still further features in the described preferred embodiments the 
method further comprises expanding a range of the second localized histogram. 

According to still further features in the described preferred embodiments the 
method further comprises linearly spreading intensity values over respective intensity 
ranges of the first and the second localized histogram. 

According to still further features in the described preferred embodiments the 
method further comprises quadratically spreading intensity values over an intensity 
range of the third localized histogram. 

According to still further features in the described preferred embodiments the 
method further comprises applying a low-pass filter on at least a portion of the 
intensity histogram, subsequently to the step of combining the plurality of improved 
localized intensity histograms. 

According to still further features in the described preferred embodiments the 
removing the clutter from the image comprises: calculating a statistical deviation for 
each picture-element over the set of still-images, thereby providing a statistical 
deviation matrix having a plurality of matrix-elements; and determining, for each 
picture element, whether a respective matrix-element of the average intensity matrix is 
above a first intensity threshold and whether a respective matrix-element of the 
statistical deviation matrix is below an additional intensity threshold, and if so then 
marking the picture- element as a clutter picture-element in the image. 



6 

According to still further features in the described preferred embodiments the 
method further comprises calculating the first intensity threshold, wherein the first 
intensity threshold is defined as an intersection point between two localized functions 
of the plurality of localized functions. 

According to still further features in the described preferred embodiments the 
method further comprises: using the statistical deviation matrix to construct a second 
intensity histogram; fitting the second intensity histogram to a second sum of a 
plurality of localized functions; and calculating the additional intensity threshold, 
wherein the additional intensity threshold is defined as an intersection point between 
two localized functions of the plurality of localized functions of the second sum. 

According to still further features in the described preferred embodiments the 
method further comprises outlining at least one region in the image, and assigning to 
each clutter picture-element an intensity value corresponding to a location of the 
clutter picture-element relative to the at least one region. 

According to still further features in the described preferred embodiments the 
outlining at least one region comprises: applying a thresholding procedure to the set of 
still-images in a Boolean manner, so as to construct at least one binary matrix having a 
plurality of binary-valued matrix -elements; and for each binary matrix of the at least 
one binary matrix, clustering the binary matrix, so as to obtain at least one cluster of 
matrix-elements having a predetermined polarity, and marking picture-elements 
corresponding to at least a portion of matrix-elements enveloping the at least one 
cluster as outlining picture-elements; thereby outlining the at least one region. 

According to another aspect of the present invention there is provided a 
method of detecting clutter in a set of images arranged grid-wise in a plurality of 
picture-elements, each image being represented by a plurality of intensity values over 
the grid, the method comprising: (a) calculating a set-averaged intensity value for each 
picture-element, thereby providing an average intensity matrix having a plurality of 
matrix-elements; (b) calculating a statistical deviation for each picture-element over 
the set of images, thereby providing a statistical deviation matrix having a plurality of 
matrix-elements; (c) determining, for each picture element, whether a respective 
matrix-element of the average intensity matrix is above a first intensity threshold and 
whether a respective matrix-element of the statistical deviation matrix is below an 
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additional intensity threshold, and if so then marking that the picture-element as a 
clutter picture-element in the set of images. 

According to further features in preferred embodiments of the invention 
described below, the method further comprises, prior to step (c): constructing a first 
5 intensity histogram characterizing the average intensity matrix, and constructing a 
second intensity histogram characterizing the statistical deviation matrix. 

According to still further features in the described preferred embodiments the 
method further comprises, prior to step (c): fitting the first intensity histogram to a first 
sum of a plurality of localized functions and using the plurality of localized functions 
10 of the first sum for calculating the first intensity threshold; and fitting the second 
intensity histogram to a second sum of a plurality of localized functions and using the 
plurality of localized functions of the second sum for calculating the second intensity 
threshold. 

According to still further features in the described preferred embodiments the 
15 plurality of localized functions of the first sum comprises a first localized function, a 

second localized function and a third localized function. 

According to still further features in the described preferred embodiments the 

first intensity threshold equals an intersection point between the second and the third 

localized functions of the first sum. 
20 According to still further features in the described preferred embodiments the 

plurality of localized functions of the second sum comprises a first localized function, 

a second localized function and a third localized function. 

According to still further features in the described preferred embodiments the 

additional intensity threshold equals an intersection point between the first and the 
25 second localized functions of the second sum. 

According to yet another aspect of the present invention there is provided a 

method of outlining at least one region in a set of images arranged grid- wise in a 

plurality of picture-elements, each image being represented by a plurality of intensity 

values over the grid and characterized by an intensity histogram, the method 
30 comprising: (a) calculating a set-averaged intensity value for each picture-element, 

thereby providing an average intensity matrix having a plurality of matrix-elements, 

and constructing a first intensity histogram characterizing the average intensity matrix; 

(b) fitting the first intensity histogram to a sum of a plurality of localized functions, so 
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as to provide at least one intensity threshold, each intensity threshold of the at least 
one intensity threshold being defined as an intersection point between two localized 
functions of the plurality of localized functions; (c) for each intensity threshold of the 
at least one intensity threshold, applying a thresholding procedure to the set of images 
5 in a Boolean manner, so as to construct at least one binary matrix having a plurality of 
binary-valued matrix-elements; and (d) for each binary matrix of the at least one 
binary matrix, clustering the binary matrix, so as to obtain at least one cluster of 
matrix-elements having a predetermined polarity, and marking picture-elements 
corresponding to at least a portion of matrix-elements enveloping the at least one 
10 cluster as outlining picture-elements; thereby outlining the at least one region. 

According to further features in preferred embodiments of the invention 
described below, the method further comprises removing clutter from at least one 
image of the set of images. 

According to still further features in the described preferred embodiments the 
15 removing the clutter from the image comprises: calculating a statistical deviation for 
each picture-element over the set of still-images, thereby providing a statistical 
deviation matrix having a plurality of matrix-elements and constructing a second 
intensity histogram characterizing the statistical deviation matrix; fitting the second 
intensity histogram to a second sum of a plurality of localized functions, so as to 
20 provide at least one additional intensity threshold, the at least one additional intensity 
threshold being defined as an intersection point between two localized functions of the 
second sum; and determining, for each picture element, whether a respective matrix- 
element of the average intensity matrix is above the second intensity threshold and 
whether a respective matrix-element of the statistical deviation matrix is below one of 
25 the at least one additional intensity threshold, and if so then marking the picture- 
element as a clutter picture-element in the image. 

According to still further features in the described preferred embodiments the 
calculating the statistical deviation for each picture-element over the set of images, 
comprises calculating a mean-square-error of a respective matrix element of the 
30 average intensity matrix. 

According to still further features in the described preferred embodiments the 
method further comprises performing at least one morphological operation on the at 
least one binary matrix. 



According to still further features in the described preferred embodiments the 
method further comprises, for each region of the at least one region, defining an origin 
of the grid, the origin being defined as a central picture-element of the region. 

According to still further features in the described preferred embodiments the 
5 method further comprises transforming the grid into a polar representation, using the 
origin of the grid, the polar representation being represented by a radial matrix and a 
angular matrix. 

According to still further features in the described preferred embodiments the 
method further comprises digitizing the radial matrix and the angular matrix using a 

10 predetermined resolution. 

According to still further features in the described preferred embodiments the 
at least one intensity threshold comprises a first intensity threshold and a second 
intensity threshold, hence the at least one binary matrix comprises a first binary matrix 
and a second binary matrix. 

15 According to still further features in the described preferred embodiments 

enveloping matrix elements of the first binary matrix corresponds to an inner boundary 
of the left ventricle, and enveloping matrix elements of the second binary matrix 
corresponds to an outer boundary of the left ventricle. 

According to still further features in the described preferred embodiments the 

20 method further comprises, for each binary matrix of the at least one binary matrix: 
identifying matrix-elements representing the papillary muscles and removing the 
identified matrix-elements from the binary matrix. 

According to still further features in the described preferred embodiments the 
identification of matrix-elements representing the papillary muscles is effected by a 

25 region-growing procedure. 

According to still further features in the described preferred embodiments the 
method further comprises for each binary matrix of the at least one binary matrix: 
identifying matrix-elements representing cusps of the mitral valve, and removing the 
identified matrix-elements from the binary matrix. 

30 According to still further features in the described preferred embodiments the 

identification of the matrix-elements representing the cusps of the mitral valve 
comprises: calculating, for each picture-element of the plurality of picture-elements of 
the grid, a deviation of an intensity value of the picture-element; and identifying 



10 

matrix-elements corresponding to picture-elements having a high deviation of 
intensity value as matrix-elements representing the cusps of the mitral valve. 

According to still further features in the described preferred embodiments the 
method further comprises: applying a set of temporal and spatial filters, so as to reject 
5 at least a few matrix-elements enveloping the at least one cluster of matrix-elements; 
and assigning to picture-elements corresponding to the rejected matrix-elements 
interpolated intensity values. 

According to still another aspect of the present invention there is provided an 
apparatus for improving an image by transforming an intensity histogram thereof, the 

10 apparatus comprising: a fitter, for fitting the intensity histogram to a sum of a plurality 
of localized functions; a histogram definer, for defining a plurality of localized 
intensity histograms using the plurality of localized functions; and a histogram 
transformer, supplemented by an algorithm for performing at least one image 
enhancement procedure, for enhancing each localized intensity histogram, thereby to 

15 provide a plurality of improved localized intensity histograms, and combining the 
plurality of improved localized intensity histograms. 

According to further features in preferred embodiments of the invention 
described below, each localized intensity histogram of the plurality of localized 
intensity histograms is characterized by an intensity range having a minimal intensity 

20 value and a maximal intensity value, such that at least one of the minimal and maximal 
intensity values coincides with an intersection point between two localized function of 
the plurality of localized functions. 

According to still further features in the described preferred embodiments the 
at least one image enhancement procedure is selected so as to enlarge a relative 

25 portion of a high-intensity region of the intensity histogram. 

According to still further features in the described preferred embodiments each 
image enhancement procedure of the at least one image enhancement procedure is 
independently selected from the group consisting of histogram equalization and 
histogram specification. 

30 According to still further features in the described preferred embodiments the 

plurality of localized functions comprises a first localized function, a second localized 
function and a third localized function, and further wherein the plurality of localized 
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histograms comprises a first localized histogram, a second localized histogram and a 
third localized histogram. 

According to still further features in the described preferred embodiments the 
histogram equalization is performed on the first localized histogram and the second 
5 localized histogram, and the histogram specification is performed on the third 
localized histogram. 

According to still further features in the described preferred embodiments the 
transformer is operable to reduce a range of the first localized histogram. According 
to still further features in the described preferred embodiments the transformer is 

10 operable to expand a range of the second localized histogram. According to still 
further features in the described preferred embodiments the transformer is operable to 
linearly spread intensity values over respective intensity ranges of the first and the 
second localized histogram. According to still further features in the described 
preferred embodiments the transformer is operable to quadratically spread intensity 

15 values over an intensity range of the third localized histogram. According to still 
further features in the described preferred embodiments the transformer is operable to 
apply a low-pass filter on at least a portion of the intensity histogram. According to 
still further features in the described preferred embodiments the low-pass filter is 
selected from the group consisting of a binomial filter and a Gaussian filter. 

20 According to still further features in the described preferred embodiments the 

image is a moving image characterized by a plurality of picture-elements, the moving 
image being formed of a set of still-images. 

According to still further features in the described preferred embodiments the 
apparatus further comprises a histogram constructor for constructing the intensity 

25 histogram of the image. 

According to still further features in the described preferred embodiments the 
histogram constructor comprises an average calculator, for calculating, for each 
picture-element of the plurality of picture-elements, a set-averaged intensity value, 
thereby to provide an average intensity matrix representing the moving image; and to 

30 construct the intensity histogram using the average intensity matrix. 

According to still further features in the described preferred embodiments the 
apparatus further comprises a preprocessing unit, for performing at least one 
preprocessing operation on the set of images. 
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According to still further features in the described preferred embodiments the 
preprocessing unit is operable to remove clutter from the image. 

According to still further features in the described preferred embodiments the 
preprocessing unit comprises: a statistical deviation calculator, for each picture- 
5 element over the set of still-images, thereby to provide a statistical deviation matrix 
having a plurality of matrix-elements; and electronic-calculation functionality for 
determining, for each picture element, whether a respective matrix-element of the 
average intensity matrix is above a first intensity threshold and whether a respective 
matrix-element of the statistical deviation matrix is below an additional intensity 
10 threshold, and if so then marking the picture-element as a clutter picture-element in the 
image. 

According to still further features in the described preferred embodiments the 
preprocessing unit further comprises an intensity threshold calculator, for calculating 
the first intensity threshold, wherein the first intensity threshold is defined as an 
15 intersection point between two localized functions of the plurality of localized 
functions. 

According to still further features in the described preferred embodiments the 
histogram constructor is operable to construct a second intensity histogram using the 
statistical deviation matrix. 
20 According to still further features in the described preferred embodiments the 

fitter is operable to fit the second intensity histogram to a second sum of a plurality of 
localized functions. 

According to still further features in the described preferred embodiments the 
intensity threshold calculator is operable to calculate the additional intensity threshold, 
25 wherein the additional intensity threshold is defined as an intersection point between 
two localized functions of the plurality of localized functions of the second sum. 

According to still further features in the described preferred embodiments the 
preprocessing unit further comprises: an outliner, for outlining at least one region in 
the image; and intensity value assigner, for assigning to each clutter picture-element an 
30 intensity value corresponding to a location of the clutter picture-element relative to the 
at least one region. 

According to still further features in the described preferred embodiments the 
outliner comprises: a thresholding unit, for applying a thresholding procedure to the 
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set of still-images in a Boolean manner, such that at least one binary matrix having a 
plurality of binary-valued matrix-elements is constructed; and a clustering unit, for 
clustering each binary matrix of the at least one binary matrix, so as to obtain at least 
one cluster of matrix-elements having a predetermined polarity, the clustering unit 
being operable to mark picture-elements corresponding to at least a portion of matrix- 
elements enveloping the at least one cluster as outlining picture-elements. 

According to an additional aspect of the present invention there is provided an 
apparatus for detecting clutter in a set of images arranged grid-wise in a plurality of 
picture-elements, each image being represented by a plurality of intensity values over 
the grid, the apparatus comprising: (a) an average calculator, for calculating a set- 
averaged intensity value for each picture-element, thereby providing an average 
intensity matrix having a plurality of matrix-elements; (b) a statistical deviation 
calculator, for calculating a statistical deviation for each picture-element over the set 
of images, thereby providing a statistical deviation matrix having a plurality of matrix- 
elements; (c) clutter identification unit, operable to determine, for each picture 
element, whether a respective matrix-element of the average intensity matrix is above 
a first intensity threshold and whether a respective matrix-element of the statistical 
deviation matrix is below an additional intensity threshold, and if so then to mark that 
the picture-element as a clutter picture-element in the set of images. 

According to further features in preferred embodiments of the invention 
described below, the statistical deviation calculator is operable to calculate a mean- 
square-error of a respective matrix element of the average intensity matrix. 

According to still further features in the described preferred embodiments the 
first and the additional intensity thresholds are predetermined. 

According to still further features in the described preferred embodiments the 
apparatus further comprises a histogram constructor, for constructing a first intensity 
histogram from the average intensity matrix and a second intensity histogram from the 
statistical deviation matrix. 

According to still further features in the described preferred embodiments the 
apparatus further comprises, a fitter, for independently fitting each of the first and 
second intensity histograms to a sum of a plurality of localized functions; and an 
intensity threshold calculator for calculating the first and the additional intensity 
thresholds, using the plurality of localized functions, wherein the first intensity 
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threshold is defined as an intersection point between two localized functions of the 
first intensity histogram, and the additional intensity threshold is defined as an 
intersection point between two localized functions of the second intensity histogram. 

According to still further features in the described preferred embodiments the 
5 plurality of localized functions of the first intensity histogram comprises a first 
localized function, a second localized function and a third localized function, whereas 
the first intensity threshold equals an intersection point between the second localized 
function and the third localized function. 

According to still further features in the described preferred embodiments the 
10 plurality of localized functions of the second intensity histogram comprises a first 
localized function, a second localized function and a third localized function, whereas 
the additional intensity threshold equals an intersection point between the first 
localized function and the second localized function. 

According to yet an additional aspect of the present invention there is provided 
15 an apparatus for outlining at least one region in a set of images arranged grid- wise in a 
plurality of picture -elements, each image being represented by a plurality of intensity 
values over the grid and characterized by an intensity histogram, the apparatus 
comprising: a histogram constructor, for constructing a first intensity histogram 
characterizing the set of images; a fitter, for fitting the first intensity histogram to a 
20 sum of a plurality of localized functions; an intensity threshold calculator, for 
calculating at least one intensity threshold, each intensity threshold of the at least one 
intensity threshold being defined as an intersection point between two localized 
functions of the plurality of localized functions; a thresholding unit, for applying a 
thresholding procedure to the set of images in a Boolean manner using each intensity 
25 threshold of the at least one intensity threshold, such that at least one binary matrix 
having a plurality of binary-valued matrix-elements is constructed; and a clustering 
unit, for clustering each binary matrix of the at least one binary matrix, so as to obtain 
at least one cluster of matrix-elements having a predetermined polarity, the clustering 
unit being operable to mark picture-elements corresponding to at least a portion of 
30 matrix-elements enveloping the at least one cluster as outlining picture-elements. 

According to still further features in the described preferred embodiments the 
apparatus further comprises an average calculator, for calculating a set-averaged 
intensity value for each picture-element, thereby to provide an average intensity matrix 
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having a plurality of matrix-elements, wherein the histogram constructor is designed to 
construct the first intensity histogram using the average intensity matrix. 

According to still further features in the described preferred embodiments the 
set of images forms a moving image. According to still further features in the 
5 described preferred embodiments the moving image comprises a cine image. 
According to still further features in the described preferred embodiments the cine 
image comprises an ultrasound cine-loop image. According to still further features in 
the described preferred embodiments the ultrasound cine-loop image comprises an 
echocardiograph cine-loop image. 

1 0 According to still further features in the described preferred embodiments each 

localized function of the plurality of localized functions is independently selected from 
the group consisting of a Gaussian function, a Lorentzian function, a hyperbolic secant 
function, a logistic distribution, a Fourier transform and a wavelet transform. 

According to still further features in the described preferred embodiments the 

15 plurality of localized functions comprises a first localized function, a second localized 
function and a third localized function. 

According to still further features in the described preferred embodiments the 
apparatus further comprises electronic-calculation functionality for performing at least 
one morphological operation on the at least one binary matrix. 

20 According to still further features in the described preferred embodiments the 

preprocessing unit is operable to remove clutter from at least one image of the set of 
images. 

According to still further features in the described preferred embodiments the 
clustering unit comprises an origin definer for defining, for each region of the at least 
25 one region, an origin of the grid, the origin being defined as a central picture-element 
of the region. 

According to still further features in the described preferred embodiments the 
clustering unit comprises a coordinate transformation unit for transforming the grid 
into a polar representation using the origin of the grid, the polar representation being 
30 represented by a radial matrix and a angular matrix. 

According to still further features in the described preferred embodiments the 
clustering unit comprises a digitizer for digitizing the radial matrix and the angular 
matrix using a predetermined resolution. 
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According to still further features in the described preferred embodiments the 
set of images forms a cine-loop image of at least the left ventricle. According to still 
further features in the described preferred embodiments the at least one region 
comprises the left ventricle. 
5 According to still further features in the described preferred embodiments the 

apparatus further comprises a preprocessing unit for performing at least one 
preprocessing operation on the set of images. 

According to still further features in the described preferred embodiments the 
preprocessing unit is operable to remove, from each binary matrix of the at least one 
10 binary matrix, matrix-elements representing the papillary muscles. According to still 
further features in the described preferred embodiments the preprocessing unit is 
supplemented by an algorithm for performing a region-growing procedure. 

According to still further features in the described preferred embodiments the 
preprocessing unit is operable to remove, from each binary matrix of the at least one 
15 binary matrix, matrix-elements representing cusps of the mitral valve. 

According to still further features in the described preferred embodiments the 
apparatus further comprises an interpolation unit for assigning interpolated intensity 
values to non-valid or removed matrix-elements, the interpolation unit being operable 
to apply a set of temporal and spatial filters prior to the assignment so as to identify 
20 the non-valid picture-elements. 

According to still further features in the described preferred embodiments the 
preprocessing unit comprises: an intensity value deviation calculator, for calculating, 
for each picture-element of the plurality of picture-elements of the grid, a deviation of 
an intensity value of the picture-element, and further wherein the matrix-elements 
25 representing the cusps of the mitral valve are correspond to picture-elements having a 
high deviation of intensity. 

The present invention successfully addresses the shortcomings of the presently 
known configurations by providing methods and apparatus for processing and 
analyzing images having performances far exceeding the prior art. 
30 Unless otherwise defined, all technical and scientific terms used herein have 

the same meaning as commonly understood by one of ordinary skill in the art to which 
this invention belongs. Although methods and materials similar or equivalent to those 
described herein can be used in the practice or testing of the present invention, suitable 
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methods and materials are described below. In case of conflict, the patent 
specification, including definitions, will control. In addition, the materials, methods, 
and examples are illustrative only and not intended to be limiting. 

Implementation of the method and system of the present invention involves 
performing or completing selected tasks or steps manually, automatically, or a 
combination thereof. Moreover, according to actual instrumentation and equipment of 
preferred embodiments of the method and system of the present invention, several 
selected steps could be implemented by hardware or by software on any operating 
system of any firmware or a combination thereof. For example, as hardware, selected 
steps of the invention could be implemented as a chip or a circuit. As software, 
selected steps of the invention could be implemented as a plurality of software 
instructions being executed by a computer using any suitable operating system. In any 
case, selected steps of the method and system of the invention could be described as 
being performed by a data processor, such as a computing platform for executing a 
plurality of instructions. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The invention is herein described, by way of example only, with reference to 
the accompanying drawings. With specific reference now to the drawings in detail, it 
is stressed that the particulars shown are by way of example and for purposes of 
illustrative discussion of the preferred embodiments of the present invention only, and 
are presented in the cause of providing what is believed to be the most useful and 
readily understood description of the principles and conceptual aspects of the 
invention. In this regard, no attempt is made to show structural details of the invention 
in more detail than is necessary for a fundamental understanding of the invention, the 
description taken with the drawings making apparent to those skilled in the art how the 
several forms of the invention may be embodied in practice. 

In the drawings: 

FIG. 1 is a schematic illustration of an apparatus for detecting clutter in a set of 
images, according to a preferred embodiment of the present invention; 

FIG. 2 is a schematic illustration of an apparatus for outlining at least one 
region in a set of images, according to a preferred embodiment of the present 
invention; 
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FIG. 3 is a schematic illustration of a clustering unit, according to a preferred 
embodiment of the present invention; 

FIG. 4 is a schematic illustration of an apparatus for improving an image or a 
set of images, according to a preferred embodiment of the present invention; 

FIG. 5 shows a three Gaussians fit, performed according to a preferred 
embodiment of the present invention to a cine-loop containing a single cardiac cycle at 
a heart rate of 65 beats/min, referred to herein after as case 1.1; 

FIG. 6 shows a division of the histogram of the first frame of case 1.1, 
according to a preferred embodiment of the present invention; 

FIG. 7 is a flowchart diagram of a method of improving an image, according to 
a preferred embodiment of the present invention; 

FIG. 8 shows a look up table calculated for case 1.1, according to a preferred 
embodiment of the present invention; 

FIG. 9 shows the result of applying the look up table of Figure 8 to the gray- 
level histogram and to each of its three Gaussian components, according to a preferred 
embodiment of the present invention; 

FIGs. lOa-f show ultrasound images of case 1.1, where Figures 10a and lOd 
show the original end-diastolic (ED) and end-systolic (ES) frames, respectively, 
Figures 10b and lOe show, respectively, the ED and ES frames after applying 
conventional histogram equalization and rejection at threshold 38; and Figures 10c and 
1 Of show, respectively, the ED and ES frames after applying the method according to 
a preferred embodiment of the present invention; 

FIG. 11 is a graph demonstrating the processing of case 1.1; 
FIG. 12 is a comparison between prior art look up tables and look up tables 
calculated for case 1 . 1 according to a preferred embodiment of the present invention; 

FIGs. 13a-f show ultrasound images of cine-loop containing a single cardiac 
cycle at a heart rate of 71 beats/min referred to herein as case 1.15, where Figures 13a 
and 13d show the original ED and ES frames, respectively; Figures 13b and 13e show, 
respectively, the ED and ES frames after applying conventional histogram equalization 
and rejection at threshold 36; and Figures 13c and 13f show, respectively, the ED and 
ES frames after applying the method according to a preferred embodiment of the 
present invention; 
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FIG. 14 is a graph demonstrating the processing of case 1.15, according to a 
preferred embodiment of the present invention; 

FIG. 15 is a comparison between prior art look up tables and look up tables 
calculated for case 1.15 according to a preferred embodiment of the present invention; 
5 FIGs. 16a-f show ultrasound images of a cine-loop containing a single cardiac 

cycle at a heart rate of 67 beats/min referred to herein as case 1.21, where Figures 1 6a 
and 16d show the original ED and ES frames, respectively; Figures 16b and 16e show, 
respectively, the ED and ES frames after applying conventional histogram equalization 
and rejection at threshold 36; and Figures 16c and 16f show, respectively, the ED and 
10 ES frames after applying the method according to a preferred embodiment of the 
present invention; 

FIG. 17 is a graph demonstrating the processing of case 1.21, according to a 
preferred embodiment of the present invention; 

FIG. 18 is a comparison between prior art look up tables and look up tables 
15 calculated for case 1.21 according to a preferred embodiment of the present invention; 

FIG. 19 is a flowchart diagram of a method of detecting clutter, according to a 
preferred embodiment of the present invention; 

FIG. 20 shows an image of an average density matrix for case 1.1., 
interchangeably referred to herein as case 2.1 ; 
20 FIGs. 21a-f show the first three frames of case 2.1 before (Figures 21a-c) and 

after (Figures 21d-f) subtraction of the average intensity matrix from each frame; 

FIG. 22 shows an ES frame for case 2.1, according to a preferred embodiment 
of the present invention; 

FIG. 23 shows an image of the average intensity matrix for case 2.1, obtained 
25 according to a preferred embodiment of the present invention; 

FIG. 24 shows an image of a statistical deviation matrix for case 2.1, obtained 
according to a preferred embodiment of the present invention; 

FIG. 25 shows a grey-level histogram obtained from the average intensity 
matrix of case 2. 1 , according to a preferred embodiment of the present invention; 
30 FIG. 26 shows a grey-level histogram obtained from the statistical deviation 

matrix of case 2.1, according to a preferred embodiment of the present invention; 
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FIG. 27 shows an image obtained by a pixel by pixel multiplication of the 
average intensity matrix by a clutter binary matrix for case 2.1, according to a 
preferred embodiment of the present invention; 

FIG. 28 shows the ES frame for case 1.15 (interchangeably referred to herein 
as case 2.2), according to a preferred embodiment of the present invention; 

FIG. 29 shows an image of the average intensity matrix for case 2.2, obtained 
according to a preferred embodiment of the present invention; 

FIG. 30 shows an image of statistical deviation matrix for case 2.2, obtained 
according to a preferred embodiment of the present invention; 

FIG. 31 shows the grey-level histogram obtained from the average intensity 
matrix of case 2.2, according to a preferred embodiment of the present invention; 

FIG. 32 shows the grey-level histogram obtained from the statistical deviation 
matrix of case 2.2, according to a preferred embodiment of the present invention; 

FIG. 33 shows an image obtained by a pixel by pixel multiplication of the 
average intensity matrix by the clutter binary matrix for case 2.2, according to a 
preferred embodiment of the present invention; 

FIG. 34 an ES frame for a cine-loop containing a single cardiac cycle at a heart 
rate of 74 beats/min referred to herein as case 2.13, according to a preferred 
embodiment of the present invention; 

FIG. 35 shows an image of the average intensity matrix for 2.13, obtained 
according to a preferred embodiment of the present invention; 

FIG. 36 shows an image obtained by a pixel by pixel multiplication of the 
average intensity matrix by the clutter binary matrix for case 2.13, according to a 
preferred embodiment of the present invention; 

FIG. 37 is a flowchart diagram of a method of outlining a region, according to 
a preferred embodiment of the present invention; 

FIGs. 38a-b are representative examples of images of two binary matrices, 
according to a preferred embodiment of the present invention; 

FIG. 39 is a representative example of an image of a binary matrix defining the 
location of the papillary muscles, according to a preferred embodiment of the present 
invention; 
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FIGs. 40a-b show an image of the average intensity matrix (Figure 40a) and 
the mitral valve region-of-interest (Figure 40b), according to a preferred embodiment 
of the present invention; 

FIGs. 41a-b are schematic illustrations of a preferred process for determining 
the left ventricle boundaries, according to a preferred embodiment of the present 
invention; 

FIGs. 42a-e show inner and outer outlines of the left ventricle performed 
according to a preferred embodiment of the present invention on case 1.1, 
interchangeably referred to herein as case 3.1; and 

FIGs. 43a-e show inner and outer outlines of the left ventricle performed 
according to a preferred embodiment of the present invention on case 1.15, 
interchangeably referred to herein as case 3.2. 

DESCRIPTION OF THE PREFERRED EMBODIMENTS 

The present invention is of methods and apparatus which can be used for 
analyzing and processing images. Specifically, the present invention can be used to 
analyze and improve the quality of still as well as moving images, including, without 
limitation, single-frame or cine-loop ultrasound images, obtained, e.g., during 
echocardiography. 

The principles and operation of methods and apparatus according to the present 
invention may be better understood with reference to the drawings and accompanying 
descriptions. 

Before explaining at least one embodiment of the invention in detail, it is to be 
understood that the invention is not limited in its application to the details of 
construction and the arrangement of the components set forth in the following 
description or illustrated in the drawings. The invention is capable of other 
embodiments or of being practiced or carried out in various ways. Also, it is to be 
understood that the phraseology and terminology employed herein is for the purpose 
of description and should not be regarded as limiting. 

According to one aspect of the present invention there is provided an apparatus 
10 for detecting clutter in a set of images. The set of images can be, for example, an 
ordered set (a series) of images, in which case the set forms a moving image, such as, 
but not limited to, a cine image. More preferably the set of images forms an 
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ultrasound cine-loop image, e.g., an echo cardiograph cine-loop image. The set is 
preferably arranged grid-wise in a plurality of picture-elements, e.g., pixels, 
arrangements of pixels, and the like. 

The term "pixel" is sometimes abbreviated herein to indicate a picture element. 
5 However, this is not intended to limit the meaning of the term "picture element" which 
refers to a unit of the composition of an image. 

According to a preferred embodiment of the present invention the grid is a 
rectangular {e.g., Cartesian) grid, but this need not necessarily be the case, as other 
grids (triangular, polar, etc.) can be defined. The set of images is preferably 
10 represented by a plurality of intensity values {e.g. , grey-levels) over the grid. 

Referring now to the drawings, Figure 1 is a schematic illustration of apparatus 
10, which, according to various exemplary embodiments of the present invention, 
comprises an average calculator 12, a statistical deviation calculator 14 and a clutter 
identification unit 16, which is preferably in communication with calculators 12 and 
15 14. Average calculator 12 calculates a set-averaged intensity value for each picture- 
element of the grid. The set-averaged intensity values are preferably represented as an 
intensity matrix, M, preferably defined as: 

M(m,n)=^Y.F p {m,n) (EQ. 1) 

p=\ 

where, P denotes the size of the set, F P denotes the /?th image of the set, and the pair 
20 {m,n) identifies the picture-element over the grid. Preferably {m,n) denotes the address 

of the picture-element (for a Cartesian grid, for example, {m,n) can represent the x-y 

location of the picture-element). 

Statistical deviation calculator 14 calculates a statistical deviation for each 

picture-element over the set. According to a preferred embodiment of the present 
25 invention calculator 14 communicates with calculator 12. Similarly to the set-average 

above, the statistical deviations are preferably represented as an intensity matrix, O. 

According to a preferred embodiment of the present invention the statistical 

deviations of the picture elements are proportional to their mean-square- errors with 

respect to the matrix-elements of M. Specifically, denoting the coefficient of 
30 proportionality by a: 
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0(m i n) = aj±f\F p (m it i)-M{m > n)] 2 . (EQ. 2) 

Preferably, the value of a is selected such that each matrix-element of O is below the 
maximal allowed intensity value. For example, for an 8-bit grey-scale image, a can be 
selected such that 0(m,n) < 255 for all matrix-elements of O. 

Clutter identification unit 16 preferably determines, for each picture element, 
hence also for each matrix-element of the matrices M and O, whether M(m,n) is above 
a first intensity threshold, T M , and whether 0(m,n) is below an additional intensity 
threshold, 7b According to a preferred embodiment of the present invention, (m,n) is 
marked as clutter picture-element if both conditions are fulfilled, i.e., M(m,n) > T M and 
0(m,n) < T 0 . 

The intensity thresholds T u and 7b can be either predetermined (e.g., supplied 
by the user) or, more preferably, can be automatically determined by apparatus 10. 
Hence, while conceiving the present invention it has been hypothesized and while 
reducing the present invention to practice it has been realized that T M and 7b can be 
found using the intensity histogram of the image. Specifically, the values of T M and 
T 0 can be found by determining the functional dependence of the intensity histogram, 
and extracting the values of T M and 7b therefrom. 

Thus, according to a preferred embodiment of the present invention apparatus 
10 comprises a histogram constructor 17, for constructing a first intensity histogram, 
H M , from the average intensity matrix, M, and a second intensity histogram, H 0 , from 
the statistical deviation matrix, O. The intensity histograms can be in any form known 
in the art such as a graph showing the number of matrix-elements of the respective 
matrix at each different intensity value found in therein. Alternatively, the histograms 
may be in a numeric form, e.g., two-dimensional arrays in which pairs of numbers 
represent intensity values and their corresponding occupation numbers in the matrices. 

In any event, once the intensity histograms are constructed, a mathematical 
procedure of fitting can be employed so as to determine their functional dependence 
from which T M and T 0 , can be extracted. According to various exemplary 
embodiments of the present invention, apparatus 10 comprises a fitter 18 and an 
intensity threshold calculator 19, being in communication thereamongst. Fitter 18 
preferably fits each intensity histogram to a sum (e.g., a weighted sum) of a plurality 
of localized functions, e.g., three localized functions. 
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As used herein "localized function" refers to any function having a local 
support and which is significantly suppressed far (say, at a distance of about 10 widths 
of the local support) from the local support. Representative examples of localized 
functions include, without limitation, a Gaussian function, a Lorentzian function, a 
hyperbolic secant function (also known as sech), a logistic distribution and the like. 
Additionally, the localized function can be represented as a series or an integral of 
other functions. For example, in this embodiment, the localized function can be a 
Fourier transform, a wavelet transform and the like. 

Intensity threshold calculator 19 preferably uses the localized functions to 
calculate T M and 7b, and provides their value to unit 16. According to a preferred 
embodiment of the present invention T M and T 0 are defined as intersection points 
between two localized functions. For example, when H M is fitted to three localized 
functions, say, G,, G 2 and G 3 , T u can correspond to an intersection point between G 2 
and G 3 . Similarly, when H 0 is fitted to three localized functions, say, g u g2 and g 3 
15 (which are not necessarily different from G u G 2 and G 3 ), T 0 can correspond to an 
intersection point between g] and g 2 . 

Once the clutter picture-elements are identified, they can be removed from the 
image or, more preferably, be assigned with interpolated intensity values. A 
representative example of clutter removal procedure is provided hereinunder in the 
20 Examples section. 

According to another aspect of the present invention there is provided an 
apparatus 20 for outlining at least one region in a set of images. The set of images can 
be of any of the aforementioned forms. Apparatus 20 is particularly useful for 
echocardiograph cine-loop image, in which case apparatus 20 can be used for outlining 
the left ventricle. 

Outlining of the left ventricle has many advantages and potential uses. For 
example, the outline of the inner boundary of the left ventricle (the endocardial wall) 
can be used for evaluating the imaged plane area of the left ventricular cavity during 
the cardiac cycle. Such evaluation allows quantitative assessment of the left 
ventricle's volume as a function of time, hence also automatic calculation of basic 
parameters, such as ejection fraction and cardiac output [Lima et al, "Diagnostic 
Imaging in Clinical Cardiology," London: Martin Dunitz, pp. 2-6, 1998]. 
Furthermore, tracking the outline of the endocardial wall enables the estimation of 
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local wall-motion velocity as a function of time, which is an indicator for muscle 
viability. 

The outline of the outer boundary (epicardial wall) can be used for estimating 
the cardiac mass and for calculating local muscle-width per timeframe, which in turn 
can be used for detecting pathologic states such as Hypertrophy. 

Reference is now made to Figure 2, which is a schematic illustration of 
apparatus 20. Apparatus 20 preferably comprises a histogram constructor which may 
be for example, histogram constructor 17. The histogram constructor serves for 
constructing a first intensity histogram which characterizes the set of images. 
According to a preferred embodiment of the present invention the histogram 
constructor calculates H M . Apparatus 20 preferably comprises a fitter, e.g., fitter 18, 
for fitting the first intensity histogram to a sum of a plurality of localized functions, 
and an intensity threshold calculator, e.g., threshold calculator 19, for calculating at 
least one intensity threshold using the localized functions. The intensity threshold(s) 
are preferably defined as intersection points between localized functions, as further 
detailed hereinabove. Each intensity threshold is preferably used for outlining a 
different region of the set or a different portion of the same region. For example, 
when the left ventricle is outlined, three localized functions, G,, G 2 and G 3 , and two 
thresholds, T x and T 2 , are preferably used, whereby T x is used for outlining the 
endocardial wall and T 2 is used for outlining the epicardial wall. In this embodiment, 
T\ preferably corresponds to the intersection point between Gi and G 2 , and T 2 
preferably corresponds to the intersection point between G 2 and G 3 . 

According to a preferred embodiment of the present invention apparatus 20 
further comprises a thresholding unit 22, for applying a thresholding procedure to the 
set of images. The thresholding is preferably applied for each intensity threshold of at 
least one intensity threshold. Additionally, the thresholding is preferably applied in a 
Boolean manner, such that one or more binary matrices having a plurality of binary- 
valued are constructed. For example, when two thresholds, T\ and T 2 are used, two 
binary matrices, BP X and & 2 can be constructed for each image, p, as follows: 
B?{m,n) = [F p ( mi n)>T,} 

VM = \f,M>T 2 \ (EQ ' 3) 
Once constructed, the binary matrices are preferably smoothed, for example, 
by applying one or more morphological operations, such as opening and closing [to 
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this effect see R. C. Gonzalez and R. W. Woods, "Digital Image Processing", pp. 518- 
548, Addison- Wesley Publishing Company, 1992]. Broadly speaking, morphological 
operation involved the passing of a convolution kernel, also known as the structuring 
element, over the binary matrices, and, based upon the convolution between the 
structuring element and the binary matrix, determining whether or not to change the 
polarity of a particular matrix-element. The morphological opening operation 
generally eliminates thin protrusions, smoothes outward bumps and breaks narrow 
sections of the matrix. The morphological closing operation generally eliminates 
small holes and removes inward bumps. 

In various exemplary embodiments of the present invention apparatus 20 
further comprises a clustering unit 24. Clustering unit 24 clusters each binary matrix 
to obtain at least one cluster of matrix-elements of a predetermined polarity. For 
example, if the binary-valued matrix-elements of the binary matrices can assume O's 
and l's, clustering unit 24 can obtain a cluster of, say, l's. According to a preferred 
embodiment of the present invention picture-elements of the set which correspond to 
matrix-elements enveloping the cluster(s) are marked as outlining picture-elements. 
The outlining picture-elements can thereafter be emphasized on each image of the set, 
e.g., using a different color or grey-level so as to distinguish them from other picture- 
elements. 

Reference is now made to Figure 3, which is a schematic illustration of 
clustering unit 24, according to a preferred embodiment of the present invention. 
Hence, clustering unit 24 preferably comprises an origin definer 32 for defining a 
central picture-element of a region, interchangeably referred to hereinafter as origin. 
The definition of the origin can also be automatic or can be inputted from an external 
source, e.g., supplied by the user. Additionally, the definition of the origin can be a 
combination of inputting and automatic definition, whereby the input is treated as a 
preliminary origin which is used for calculating the central picture-element. 
According to a preferred embodiment of the present invention, the definition of the 
origin is performed once on the average intensity matrix, M, and subsequently used for 
all the images in the set. 

Clustering unit 24 preferably comprises a coordinate transformation unit 34 for 
transforming the grid of the images into a polar representation using the origin. As 
will be appreciated by one ordinarily skilled in the art, the polar representation of each 
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image, intensity matrix or binary matrix can comprise a radial matrix, R, storing the 
distance of a respective picture- or matrix-element from the origin, and an angular 
matrix, 6, storing the angle of a respective picture- or matrix-element relative to the 
one of the axes of the original coordinate system. According to a preferred 
embodiment of the present invention clustering unit 24 further comprises a digitizer 36 
which digitizes the matrices R and 0, in accordance with the to the required axial and 
lateral resolution. 

The use of polar representation facilitates the clustering of the images. 
Specifically, the clustering can be done by traversing the matrices on a radial direction, 
e.g., from the origin outwards, and determining, for each angle stored in 0, when the 
polarity of the matrix element is inverted. For example, supposing that the matrix- 
elements corresponding to the background of the image are assigned with zeroes, 
while non-background matrix-elements are assigned with non-zeroes. Then, for a 
given angle, clustering unit 24 traverses the radial direction, locates the first non-zero 
matrix element and marks the corresponding picture-element as an outlining picture- 
element of the inner boundary of the region. Continuing the traverse along the radial 
direction, clustering unit 24 preferably locate the next zero element and marks the 
corresponding picture-element as an outlining picture-element. This procedure can be 
repeated for all angles stored in 6> such that the inner and/or outer boundary of the 
region is outlined. 

The outlining of the region can be preceded by one or more preprocessing 
operations, for improving the efficiency of the clustering procedure, minimizing 
misclassifications and/or improving the quality of the images. Thus, according to a 
preferred embodiment of the present invention, apparatus 20 comprises a 
preprocessing unit 23, for performing at least one preprocessing operation. Many 
preprocessing operations and combination thereof are contemplated, representative 
examples include, without limitation, removal, replacement and interpolation of 
picture-elements, either temporarily, e.g., for improving the efficiency of the clustering 
process, or permanently, e.g., for improving the quality of the resulting image. 

Typically, the preprocessing is directed at handling picture-elements identified 
as storing no or irrelevant information. Such picture-elements can include, without 
limitation, clutter picture-elements and picture-element identified as not belonging to 
the outlined region. For example, when apparatus 20 is used for outlining the left 
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ventricle in an ultrasound image, picture-elements corresponding to the papillary 
muscles can obstruct the processing procedure, because the papillary muscles 
oftentimes appear on the ultrasound image as high-intensity regions within the left 
ventricular cavity, hence may lead to misclassification of the corresponding picture- 
elements. Additionally, picture-elements corresponding to the base of the mitral valve 
can contribute to discontinuity in the appearance of the boundary. 

Hence, one operation which can be performed by preprocessing unit 23 is 
clutter removal. Specifically, preprocessing unit preferably identifies and removes 
clutter picture-elements from all or a few images of the set. This operation can be 
achieved by associating preprocessing unit 23 with apparatus 10. 

Another operation which can be performed by preprocessing unit 23 is the 
removal of matrix -elements representing the papillary muscles. According to a 
preferred embodiment of the present invention the identification of papillary muscles 
is done by applying the thresholds T { and T 2 on the average intensity matrix, M, in a 
Boolean manner to construct two (average) binary matrices, B x and B 2i and 
performing a region-growing procedure on B l and B 2 . Region-growing techniques are 
known in the art and are found, e.g., in an article by S. W. Zucker, entitled "Region 
Growing: Childhood and Adolescence," published in Computer Graphics and Image 
Processing, vol. 5, pp. 382-399, 1976. A detailed example for the identification and 
removal of matrix-elements corresponding to papillary muscles is provided in the 
Examples section that follows. 

An additional operation which can be performed by preprocessing unit 23 is 
the removal of matrix-elements representing cusps of the mitral valve. The 
identification of these matrix-elements is preferably based on the intensity value 
deviation, A of the corresponding picture-elements. Specifically, matrix-elements 
representing cusps of the mitral valve preferably correspond to picture-elements 
having a high deviation of intensity. The deviation of intensity for the pth image is 
preferably calculated using the average intensity matrix, M, as follows: 

D(m ' n) M(m, n) ' ^Q- 4 ) 

The present embodiments thus successfully handle the appearance of picture- 
elements corresponding to no or irrelevant information by removing them from the 
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image, and, optionally and preferably, assigning them with, e.g., interpolated intensity 
values, prior to the clustering procedure. 

Optionally and preferably, apparatus 20 further comprises an interpolation unit 
26 for assigning interpolated intensity values to non-valid matrix-elements. Such non- 

5 valid matrix-elements can correspond, e.g., to clutter, papillary muscles, cusps of the 
mitral valve or any other type of information, for which interpolation is needed. For 
example, a set of temporal and spatial filters can be used to identify non-valid matrix- 
elements. As further demonstrated in the Examples section that follows, the 
interpolation is particularly useful for improving the calculated outlines of the region. 

10 For example, the interpolation can be used for handling non-valid or incorrect 
estimates of the outline location in different regions of different frames. 

According to an additional aspect of the present invention there is provided an 
apparatus 40, for improving an image or a set of images. Apparatus 40 can be used for 
improving one image, a set of images or a plurality of sets of images. For example, 

15 apparatus 40 can be used in analysis of echocardio graph procedure of different 
patients or different echocardio graph procedures of the same patient, whereby for each 
cine-loop, the analysis is independent of the other cine-loops. The advantage of this 
embodiment is that apparatus 40 automatically tailors a histogram transformation to 
each specific cine-loop, of each patient. 

20 Reference is now made to Figure 4 which is a schematic illustration of 

apparatus 40. Apparatus 40 preferably comprises a fitter, e.g., fitter 18, for fitting an 
intensity histogram, H, to a sum of a plurality of localized functions, as further 
detailed herein above. H can be either the intensity histogram of the image or an 
intensity histogram representing the set of images, e.g., H M , the intensity histogram of 

25 the average intensity matrix, M. Apparatus 40 further comprises a histogram definer 
42, which uses the localized functions for defining a plurality of localized intensity 
histograms, denoted H u where /is an integer running from 1 to the number of 
localized intensity histograms. Typically, but not obligatory, three localized intensity 
histograms are defined. 

30 According to a preferred embodiment of the present invention each localized 

intensity histogram is characterized by an intensity range defined using intersection 
points between two localized functions. For example, suppose that H{u) is fitted to the 
sum of three localized functions, G\, G 2 and G 3 , whereby G\ intersects with G 2 at 
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u = Ti and G 2 intersects with G 3 at u = r 2 , where F 2 > T h Suppose further than the 
intensity range of H(u) is from u = Wmin to u = u max (for 8-bit image w min can be 0 and 
w max can be 255). Then, according to a preferred embodiment of the present invention 
histogram definer 42 defines three localized intensity histograms, H u H 2 and H 3 , such 

5 that Hi spans from a minimal intensity value of w min to a maximal intensity value of T x ; 
#2 spans from a minimal intensity value of T\+\ to a maximal intensity value of T 2 ; 
and H 3 spans from a minimal intensity value of F 2 +l to a maximal intensity value of 
w max . Each localized intensity histogram can be padded by zeroes if desired. 

Apparatus 40 further comprises a histogram transformer 44 which enhances 

10 each localized intensity histogram, Hi, and combines the resultant, improved, localized 
intensity histogram into a single histogram. The enhancement of the H\ is preferably 
selected so as to enlarge a relative portion of a high-intensity region of the original 
intensity histogram, H This can be done by any technique known in the art, such as, 
but not limited to, histogram equalization, histogram specification and the like. When 

15 three localized intensity histogram are defined, H, and H 2 are preferably subjected to 
histogram equalization and H 3 is preferably subjected to histogram specification. 

Prior to the combining of the His they can be subjected to further operations 
such as reduction or expansion. More specifically, the range of H\ is preferably 
reduced by a predetermined reduction factor, 0d, while the range of H 2 is preferably 

20 expanded by a predetermined expansion factor a 2 . As will be appreciated by one of 
ordinary skill in the art, once the reduction or expansion are performed, the intensity 
ranges of the His no longer coincide with the intersection points between the localized 
functions. 

Optionally and preferably, one or more of the His can be spread over its 
25 respective range, whether or not this range was subjected to reduction or expansion. 
Preferably, for three His, a linear spreading is employed on the intensity ranges of Hi 
and H 2 while a non-linear spreading, e.g., (quadratic) spread is employed on the 
intensity ranges of H 3 . 

According to a preferred embodiment of the present invention, once the His 
30 are combined, transformer 44 applies a low-pass filter (e.g., a binomial filter or a 
Gaussian filter) on at least a portion of the resulting intensity histogram. The low-past 
filter is preferably employed near the intersecting points between the localized 
functions so as to minimize information loss. 
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When used for improving a set of images or a plurality of sets of images, 
apparatus 40 can include or be associated with apparatus 10 and/or apparatus 20, so as 
to remove clutter and/or outline a particular region in the set(s). For example, when it 
is desired to perform analysis of one or more ultrasound cine-loops of the left 
5 ventricle, apparatus 40 can comprise apparatus 10 and apparatus 20. In this 
embodiment, apparatus 40, preferably eliminates clutter, outlines the endocardial 
and/or epicardial walls, and transforms the histograms of each image so as to provide 
the physician with an improved, substantially clutter-free cine-loop, vividly 
designating the left ventricle. 

10 

It is expected that during the life of this patent many relevant imaging 
techniques will be developed and the scope of the term set of images is intended to 
include all such new technologies a priori. 

15 Additional objects, advantages and novel features of the present invention will 

become apparent to one ordinarily skilled in the art upon examination of the following 
examples, which are not intended to be limiting. Additionally, each of the various 
embodiments and aspects of the present invention as delineated hereinabove and as 
claimed in the claims section below finds experimental support in the following 

20 examples. 

EXAMPLES 

Reference is now made to the following examples, which together with the 
above descriptions illustrate the invention in a non limiting fashion. 

25 EXAMPLE I 

Adaptive Transformation of Intensity Histogram 
The teachings of the present embodiments were employed for performing 
adaptive transformation of intensity histograms. As further demonstrated below, the 
application of the teachings of the present embodiments significantly improves the 
30 quality of the Echocardiograph images. 
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Materials and Methods 



Theoretical Background 

The processing of the images was based on the fitting of the gray-level 
histogram curve to a sum of three Gaussians, G\, G 2 and G 3 , representing different 
5 regions in the image. Gi represented the left ventricular cavity (low intensity 
Gaussian), G 2 represented low-intensity pixels within the cardiac muscle (medium 
intensity Gaussian), and G 3 represented high-intensity pixels within the cardiac muscle 
(high intensity Gaussian). 



10 case 1.1, which is fitted to three Gaussians. In Figure 5, the solid line correspond to 
the original histogram, the dash-dotted line correspond to the sum of the three 
Gaussian, and the dotted lines correspond to the individual Gaussian G], G 2 and G 3 . 

The histogram was then divided into three localized histograms, each 
dominated by a different Gaussian. In order to achieve maximal separation, the 

15 division points, T\ and T 2 , were set at the intersection points between adjacent 
Gaussian functions. In the case 1.1, the division points were 38 and 92, on a scale 
between 1 and 256. 

Figure 6 shows the division of the histogram of the first frame of case 1.1. The 
original image appears on the top left. Pixels whose gray level is lower than T\ (38 in 
20 this case) appear in white on the top right. Pixels whose gray level is between T\ and 
T 2 (92) appear in white on the bottom left. Pixels whose gray level exceeds T 2 are 
shown on the bottom right. 

For the purpose of measuring the separation, the following parameters have 
been defined: 



Figure 5 shows a gray-level histogram of a cine-loop, referred to hereinafter as 



25 




(EQ. 5) 



\G x {u)du 




(EQ. 6) 



J[G 2 («)+G 3 («)]rf« 



-00 
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}[G,(«)+G 2 («)]rf U 




J[G,(«)+G 2 («)]rf« 



(EQ. 7) 




lu 



(EQ. 8) 



PJ0' = 1, 2) are detection probabilities which measure the probability for a 
pixel whose intensity exceeds T t to belong to the expected Gaussian function {Gi or G 3 
for i = 1, G3 for i = 2). (/ = 1, 2) are false-alarm probabilities which measure the 
probability for a pixel whose intensity exceeds F, to belong to the incorrect Gaussian 
function (G\ for i=l, Gi or G 2 for / = 2). Ideal separation is achieved when 
P] = P 2 d =1.0 and Pj^ = Pj = 0.0 . 

In most cases, the information in the left ventricular cavity (the first histogram 
region) has little diagnostic use, while the majority of the relevant data resides within 
the cardiac muscle (including the second and third regions). The second region 
encloses many or most of the pixels (in the example shown in Figure 6, 33.7 % of the 
pixels belong to the first section, 42.3 % belong to the second section and 24.0 % 
belong to the third section), but tends to utilize only a small portion of the gray-level 
scale. Therefore, a compression was employed to the first localized histogram, while a 
stretching was employed to the second localized histogram. In addition, histogram 
specification [Jain AK, "Fundamentals of Digital Image Processing", New Jersey: 
Prentice Hall, pp. 241-244, 1989] was applied to each localized histogram separately, 
so as to improve image quality. 

Image Processing 

While examining the raw data, it was observed that the calculated division 
points change only slightly between adjacent frames of the cine-loop. Therefore, to 
this end the calculation was performed on a single histogram for the entire cine-loop. 
This calculation has the advantage of decreasing the probability for receiving singular 
gray-level distributions. 



34 



Reference is now made to Figure 7, which is a flowchart diagram of a method 
of improving an image, according to a preferred embodiment of the present invention. 
A description of the specific steps employed in the present examples accompanies the 
general description of the method. 

It is to be understood that, unless otherwise defined, the method steps 
described hereinbelow can be executed either contemporaneously or sequentially in 
many combinations or orders of execution. Specifically, the ordering of the flowchart 
of Figure 7 is not to be considered as limiting. For example, two or more method 
steps, appearing in the following description or in the flowchart of Figure 7 in a 
particular order, can be executed in a different order (e.g., a reverse order) or 
substantially contemporaneously. 

Referring now to Figure 7, the method preferably begins at Step 70 in which a 
set-average intensity value is calculated for each picture-element (a pixel in the present 
example) to provide the average intensity matrix, M (see Equation 1). According to a 
preferred embodiment of the present invention the averaging is performed to at least a 
single cardiac cycle. In Step 71, //(«), the intensity histogram of M, is preferably 
constructed. In Step 72, H is preferably fitted to a sum of localized function. In the 
present example, the following localized functions were used: 



In Step 73 the intersection points, T\ and T 2 between the localized functions are 
preferably calculated, and in Step 74 localized histograms are defined. In the present 
example, the following localized histograms were defined: 




(EQ. 9) 



y = G l {x)+G 2 (x)+G 3 (x) 
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».(«) 



H(u) 



u<T x 



0 



else 



H 2 (u) 



H(u) 



T 2 >u>T i 



(EQ. 10) 



0 



else 



H,{u) 



H(u) 



else 



Once the localized histograms are defined, the method proceeds to Step 75 in 
which the localized histograms are subjected to an image enhancement procedure. In 
the present example, H\(u) and H 2 (u) underwent histogram equalization, in a manner 
such that each localized histogram uniformly utilized the entire gray-level range. The 
histogram equalization of H\{u) and H 2 (u) was represented as look-up tables (LUTs), 
respectively denoted by L\(u) and L 2 (u). Further in the present example, H 3 (u) 
underwent histogram specification, described as a LUT denoted by L 3 (u). 

The histogram specification was designed to produce a histogram which 
substantially matches the following function: 



where p\, fh and fh are constants, and fix) was normalized to a unit area. As can be 
understood from Equation 11, /(x) is a constant up-to fh and matches a right half of a 
Gaussian thereafter. The reason of using this process on H 3 (u) (rather that histogram 
equalization), was to reduce the weight of high- intensity pixels, whose presence may 
cause the image to seem saturated. 

In Step 76, one or more of the localized histograms is preferably subjected to a 
range reduction or expansion. In the present example, the range of H x was reduced by 
a factor ct\ and H 2 was expanded by a factor a 2 . It was found that the selection 
on =2.0 and #2=1.75 was sufficient to provide satisfactory results. Once 
reduction/expansion operations were performed, the divisional points between H\, H 2 
and H 3 were: 



where C\ denotes the divisional point between H\ and H 2i and C 2 denotes the 




(EQ. 11) 




(EQ. 12) 
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divisional point between H 2 and // 3 . In Equation 12, the symbol L.J denotes the 
integer function "floor". 

In Step 77, the localized histograms are redistributed using a predetermined 
distribution procedure. In the present example, the redistribution is done using a set of 

5 LUTs as follows: a first LUT, Lc\{u\ linearly spreads the gray-levels between a 
minimal value (1) and C\, a second LUT, Lc 2 {u\ linearly spreads the gray-levels 
between {C\ + 1) and C 2 , and a third LUT, Ic 3 (w), quadratically spreads the gray- 
levels between (C 2 + 1) and the maximal value (256). The quadratic distribution was 
characterized by a decreasing derivative with respect to the intensity. The composite 

10 LUT for each localized histogram was then calculated as follows: 

L 2c2 (u)=L c2 (L 2 (u)) (EQ.13) 
L 2c2 (u) = L c3 (l 3 (u)) 

In Step 78 the localized histograms are combined to a single, improved, 
histogram. Specifically to the present example, an improved LUT, denoted L{u), was 
obtained by stitching the three composite LUTs together: 

Xi(«) u<T, 

15 L(u)= « L 2c2 (u) T 2 >u>T } (EQ. 14) 

L^{u) u>T 2 

In Step 79, a low-pass filter is applied for at least a portion of the intensity 
range so as to produce a continuous transformation to the histogram. In the present 
example, a binomial low-pass filter was used in intensity intervals which are close to 
intersecting points T\ and T 2 \ 

20 *Mr,- r> 2;+rD=AGfe[7;- r ,7;+r]> 

4u e [t 2 - r ,T 2 +r D=LA«z [t 2 -y,T 2 + rl (EQ - 15) 

where, 

4(«)=Z(w)*£[l 6 15 20 15 6 l], (EQ. 16) 

and y is a constant, selected to be 4. 

Figure 8 shows L(u\ calculated for Case 1.1. The values of T\ and T 2 are 
25 marked by dotted lines. 

Figure 9 shows the result of applying L(u) to the gray-level histogram and to 
each of its Gaussian components, G u G 2 and G 3 . Comparison between Figure 5 (the 
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original gray-level histogram of case 1.1) and Figure 9 visually demonstrates the 
processing required for the algorithm. The first section has been compressed, while 
the second section has been stretched. Additionally, as a result of applying histogram 
equalization/specification to each section, the three components of the new gray-level 
5 histogram are asymmetric functions, as opposed to the Gaussian functions in the 
original case. Since the compression and stretching are applied to sections of the 
histogram rather than to the Gaussian functions themselves, these functions also 
change their form at C\ and 

Note that the abovementioned similarity in the gray-level histograms of 
10 different frames within the cine-loop can be used to deal with long cine-loops: the 
method of the present embodiments can be calculated for the first cardiac cycle, and 
subsequently applied on-line to other cardiac cycles. 

Data acquisition 

23 randomly selected cine-loops of echocardio graph images of 10 different 
15 patients were used. All cine-loops were in apical four-chamber and apical two- 
chamber views. The data has been collected using a Vivid™ 3 imaging platform, by 
GE Medical Systems Ultrasound. 

For each loop, both the teachings of the present embodiments and conventional 
histogram equalization combined with rejection were employed. Final parameter 
20 optimization was based on expert opinion of renowned cardiologists and 
echocardiography specialists. 



Results 

Table 1 below presents the values of the detection probabilities (P d \P]) and 
25 the false-alarm probabilities (P*,Pf,P} a and Pj a \ for the 23 cine-loops, referred to as 
cases 1.1-1.23. 



Table 1 



Case 
No. 




P l 


p2 


P 2 


1.1 


0.823 


0.076 


0.689 


0.039 


1.2 


0.891 


0.054 


0.734 


0.033 


1.3 


0.894 


0.055 


0.733 


0.032 


1.4 


0.853 


0.047 


0.666 


0.033 


1.5 


0.858 


0.047 


0.664 


0.030 


1.6 


0.882 


0.050 


0.721 


0.036 
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1.7 


0.887 


0.053 


0.719 


0.035 


1.8 


0.883 


0.051 


0.722 


0.036 


: i.9 


0.887 


0.053 


0.719 


0.035 


1.10 


0.898 


0.297 


0.726 


0.054 


1.11 


0.845 


0.154 


0.770 


0.044 


1.12 


0.857 


0.197 


0.776 


0.051 


1.13 


0.926 


0.250 


0.812 


0.065 


1.14 


0.859 


0.194 


0.777 


0.050 


1.15 


0.856 


0.190 


0.776 


0.048 


1.16 


0.869 


0.196 


0.763 


0.021 


1.17 


0.874 


0.106 


0.681 


0.033 


1.18 


0.864 


0.051 


0.840 


0.107 


1.19 


0.874 


0.048 


0.732 


0.047 


1.20 


0.894 


0.083 


0.812 


0.079 


1.21 


0.905 


0.120 


0.847 


0.024 


1.22 


0.801 


0.237 


0.623 


0.037 


1.23 


0.858 


0.144 


0.722 


0.034 


Mean 


0.871 


0.119 


0.740 


0.044 


Std 


0.027 


0.079 


0.056 


0.019 



As shown in Table 1, the values of the detection probabilities are large (mean 
values of 0.871 ± 0.027 for Ti and 0.74 ± 0.056 for T 2 ) while the values of the false- 
alarm probabilities are low (0.1 19 ± 0.079 and 0.044 ± 0.019) demonstrating that the 
intensity histogram was highly correlated with the sum of three Gaussians. 

Figures 10-18 show results for three representative cine-loops: case 1.1 
(Figures 10-12) case 1.15 (Figures 13-15) and case 1.21 (Figures 16-18). These cine- 
loops are characterized by high, medium and low image quality, respectively. Case 
1.1 contained a single cardiac cycle at a heart rate of 65 beats/min, case 1.15 contained 
a single cardiac cycle at a heart rate of 71 beats/min, and case 1.21 contained a single 
cardiac cycle at a heart rate of 67 beats/min. Each cine-loop was formed of 40 to 50 
frames, of which end-diastolic (ED) and the end-systolic (ES) frames are presented 
below (the ED and ES were determined based on the left ventricle's area within the 
imaging plane). 

Figures lOa-f show ultrasound images of case 1.1, where Figures 10a and lOd 
show the original ED and ES frames, respectively; Figures 10b and lOe show, 
respectively, the ED and ES frames after applying conventional histogram equalization 
and rejection at threshold 38; and Figures 10c and lOf show, respectively, the ED and 
ES frames after applying the method according to a preferred embodiment of the 
present invention. As shown in Figures 10a-f s the results obtained according to the 
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present embodiments appear crisper, and have richer gray-levels in the important 
regions, compared to both the original and the histogram equalized images. 

Figure 11 shows five curves for the processing of case 1.1: the original 
histogram is shown as a solid line, the fitted curve is shown as a dash-dotted line, and 
5 the three calculated Gaussians are shown as dotted curves. As shown in Figure 1 1, the 
fitted curve is very similar to the original histogram, and the three Gaussians are well 
separated. 

Figure 12 shows LUTs calculated for case 1.1: the LUT used for histogram 
equalization after rejection is shown as a dash-dotted line, the LUT obtained according 

10 to the present embodiments is shown as a solid line. The values of T\ and T 2 are 
marked by dotted lines. The differences between the LUTs are vivid. 

Figures 13a-f show ultrasound images of case 1.15, where Figures 13a and 13d 
show the original ED and ES frames, respectively; Figures 13b and 13e show, 
respectively, the ED and ES frames after applying conventional histogram equalization 

15 and rejection at threshold 36; and Figures 13c and 13f show, respectively, the ED and 
ES frames after applying the method according to a preferred embodiment of the 
present invention. As shown in Figures 13a-f, the results obtained according to the 
present embodiments appear sharper compared to both the original and the histogram 
equalized images. 

20 Figure 14 shows five curves for the processing of case 1.15: the original 

histogram is shown as a solid line, the fitted curve is shown as a dash-dotted line, and 
the three calculated Gaussians are shown as dotted curves. As shown in Figure 14, the 
fitted curve is very similar to the original histogram, and that the three Gaussians are 
well separated. 

25 Figure 15 shows LUTs calculated for case 1.15: the LUT used for histogram 

equalization after rejection is shown as a dash-dotted line, the LUT obtained according 
to the present embodiments is shown as a solid line. The values of T\ and Ti are 
marked by dotted lines. The differences between the LUTs are vivid. 

Figures 16a-f show ultrasound images of case 1.21, where Figures 16a and 16d 

30 show the original ED and ES frames, respectively; Figures 16b and 16e show, 
respectively, the ED and ES frames after applying conventional histogram equalization 
and rejection at threshold 36; and Figures 16c and 16f show, respectively, the ED and 
ES frames after applying the method according to a preferred embodiment of the 
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present invention. As shown in Figures 16a-f, the results obtained according to the 
present embodiments appear sharper compared to both the original and the histogram 
equalized images. 

Figure 17 shows five curves for the processing of case 1.21: the original 
5 histogram is shown as a solid line, the fitted curve is shown as a dash-dotted line, and 
the three calculated Gaussians are shown as dotted curves. As shown in Figure 17, the 
fitted curve is very similar to the original histogram, and the three Gaussians are well 
separated. 

Figure 18 shows LUTs calculated for case 1.21: the LUT used for histogram 
10 equalization after rejection is shown as a dash-dotted line, the LUT obtained according 
to the present embodiments is shown as a solid line. The values of T\ and T 2 are 
marked by dotted lines. The differences between the LUTs are vivid. 

EXAMPLE 2 

1 5 Clutter Detection 

The teachings of the present embodiments were employed for detecting clutter 
in echocardiograph images. As further demonstrated below, the application of the 
teachings of the present embodiments allows the determination of a substantial amount 
of the clutter pixels, with relatively low false detection. 

20 

Materials and Methods 
Theoretical Background 

The clutter determination procedure was based on the fact that the cardiac wall 
motion is much faster than the motion of the surrounding tissue. During a single 
25 cardiac cycle, the organs primarily contributing to the clutter can be perceived as 
standing still with respect to the transducer. Hence, according to a preferred 
embodiment of the present invention stationary or substantially stationary pixels, are 
identified as clutter pixels. 

Image Processing 

30 According to various exemplary embodiments of the present invention the 

detection of clutter can be done by executing several method steps. Reference is now 
made to Figure 19, which is a flowchart diagram of a method of detecting clutter, 
according to a preferred embodiment of the present invention. A description of the 
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specific steps employed in the present example accompanies the general description of 
the method. As stated hereinabove, the method steps described hereinbelow can be 
executed either contemporaneously or sequentially in many combinations or orders of 
execution, and the ordering of the flowchart of Figure 19 is not intended to be limiting. 

5 Hence, the method preferably begins at Step 90 in which a set-average 

intensity value is calculated for each picture-element (a pixel in the present example) 
to provide the average intensity matrix, M. The method proceeds to optional Step 91, 
in which M is subtracted from each frame to provide a clutter-cancelled frame, C p : 

C p (m,n)=F p (m,n)-M(m,n). (EQ. 17) 

10 Note that C p (m,n) can receive both positive and negative values. C p (m,ri) may 

be displayed in more than one way. In one embodiment, C p has positive entries only, 
in another embodiment, C p has negative entries only, in an additional embodiment C p 
has absolute value entries, and in a further embodiment C p has both negative and 
positive entries. 

15 Step 91 can be used to accentuate changes from one frame to the next. This is 

highly advantageous for algorithms applying local tracking of regions-of-interest 

during the cardiac cycle. 

Optionally and preferably, the method proceeds from step 91 (or from step 90 

if step 91 is not employed) to step 92 in which a statistical deviation is calculated for 
20 each picture-element, to provide the statistical deviation matrix, O (see, e.g., Equation 

2). 

For a clutter pixel (m,n\ MQfi.n) is expected to be high (only strong clutter has 
substantial effect on the image), while 0(m,ri) is expected to be low (small changes in 
time). Thus, according to a preferred embodiment of the present invention, the 

25 method proceeds to decision step 93 in which the value of each matrix-element of Mis 
compared to the threshold T M and the value of each matrix-element of O is compared 
to the threshold T 0 . If a particular matrix-element of M is above T M and a respective 
matrix-element of O is below To the method proceeds to step 94 in which the 
corresponding picture element (pixel in the present example) is marked as a clutter 

30 picture-element. 



The above steps can be executed by defining a binary matrix, C, as follows: 

C(m i r,) = [M(m,n)>T M ]n[0{m,n)<T o l (EQ. 18) 

where C receives 1 for clutter pixels, and 0 for other pixels. 

According to various exemplary embodiments of the present invention 

5 decision step 93 is preceded by step 95 in which the thresholds T M and T 0 are 
calculated, preferably by fitting the histograms H M and H 0 characterizing M and O to a 
sum of a plurality of localized functions and defining the thresholds as intersection 
points between the localized functions. 

In the present example, the three Gaussians of Example 1 (see Equation 9) 

10 were used. For H M , a lower threshold, T XM , was set to the intersection point between 
Gi and G 2 , and a higher threshold, T 2M , was set to the intersection point between G 2 
and G 3 . Similarly, for H 0 , a lower threshold, T lo , was set to the intersection point 
between Gj and G 2 , and a higher threshold, T 2 o, was set to the intersection point 
between G 2 and G 3 . The thresholds T M and T 0 were defined by setting T M = T 2M and 

15 T 0 =Tia 

For the purpose of investigating the performance of the method, the following 
ratios were defined: a ratio, denoted P c , was defined between the number of pixels 
detected as clutter and the total number of pixels in the data-region; a second ratio, 
denoted P d , was defined between the number of pixels properly detected as clutter and 
20 the total number of clutter pixels; a third ratio, denoted, P fa , was defined between the 
number of pixels falsely detected as clutter and the total number of pixels detected as 
clutter; and a fourth ratio, denoted P fat , was defined between the number of pixels 
falsely detected as clutter and the total number of pixels in the data-region. 

Note that Pd and P& can only be evaluated visually. Thus, the accuracy of Pd 
25 has been set to ± 10 %, while the accuracy of P fa has been set to ± 5 %. Also note that 
p M = p c .p /am The effectiveness of the method can be quantified by the ratios P d and 
P fat , where high values of P d and low values of P fat correspond to an effective method. 
For an ideal method P d = 1 and Pm = 0. 
Data acquisition 

30 1 6 cine-loops of echocardiograph images recorded of 1 6 different patients were 

selected. The cine-loops were in apical four-chamber and apical two-chamber views. 
The data has been collected using a Vivid™ 3 imaging platform, by GE Medical 
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Systems Ultrasound. The capability of the present embodiment is illustrated below 
using two clutter-affected loops and one clutter-free loop. 

Results 

5 Figure 20 shows an image of the average density matrix U for case 1.1 (see 

also Example 1) interchangeably referred to herein as case 2.1. As shown, the time 
averaging smears the image, but the main cardiac features are still visible. 

Figures 21a-f show the first three frames of case 2.1 before (Figures 21a-c) and 
after (Figures 21d-f) the subtraction of M from each frame (see C p Equation 17). Note 

10 that in C p the background color is gray rather than black. This results from the fact, 
that C p may receive both positive and negative values. Zero is marked using the gray 
background color, while positive values are lighter and negative values are darker. As 
shown in Figures 21a-f, C P visually highlights temporal changes. 

Table 2 below presents the values of the ratios P c , P d , P& and P fat , for each of 

15 the 1 6 cine-loops, referred to as cases 2.1-2.16. 



Table 2 



Case No. 


Pc\%\ 


Pd \%] 


Pfa I%1 


Pfat\%] \ 


2.1 


0.07 


100 


100 


0.07 


2.2 


3.79 


100 


5 


0.19 


2.3 


0.06 | 


100 


0 


0.00 


2.4 


0.39 


80 


0 


0.00 


2.5 


0.02 


100 


100 


0.02 


2.6 


0.22 


100 


100 


0.22 | 


2.7 


0.11 


50 


5 


0.01 


2.8 


0.03 


80 


25 


0.01 


2.9 


0.30 


80 


0 


0.00 1 


2.10 


0.72 


90 


0 


0.00 


2.11 


16.55 


90 


5 


0.83 


2.12 


0.35 


100 


100 


0.35 1 


2.13 


1.09 


90 


0 


0.00 


2.14 


0.00 


100 


0 


0.00 


2.15 


0.00 


100 


0 


0.00 


2.16 


0.54 


100 


100 


0.54 j 


Mean ± Std 


1.52 ±4.12 


91.3 ±13.6 


33.8 ±46.5 


0.14 ± 0.24 



As shown the values of P d are high and the values of Pf at are low, 
demonstrating the efficiency of the present embodiment of the invention. 

Figures 22-36 show results for three representative cine-loops: case 2.1 
(Figures 22-27) case 2.2 (Figures 28-33) and case 2.13 (Figures 34-36). Case 2.1 was 
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a clutter-free cine-loop, while cases 2.2 and 2.13 were clutter-affected cine-loops. All 
three loops have been gathered using high frame-rate (75 frames/s). Each loop 
contained a single cardiac cycle with a heart rate of 65 beats/min (case 2.1), 71 
beats/min (case 2.2) and 74 beats/min (case 2.13). 

5 Figure 22 show the ED frame for case 2.1. As shown in Figure 22, the ED 

frame is substantially clutter-free. Figure 23 shows an image of M for case 2.1. As 
shown in Figure 23, the averaging procedure, although smearing the image, preserves 
most of the relevant features. Figure 24 shows an image of O for case 2.1. Note that 
pixels with high intensity indicate substantial changes in the gray-level during the 

10 cardiac cycle. Figures 25-26 show the grey-level histograms H M (Figure 25) and H 0 
(Figure 26) for case 2.1. Five curves are shown in Figures 25-26, describing the 
original histogram (solid), the fitted curve (dashed), and the three calculated Gaussians 
(dotted). The intersection points between the Gaussians are 38 and 91 for H M , and 25 
and 50 for H 0 . Thus, T M =9l and T 0 = 25. As shown in Figures 25-26 the fitted 

15 curves are similar to the original histograms, and the three Gaussians are well 
separated. Figure 27 shows an image obtained by a pixel by pixel multiplication of M 
by the clutter binary matrix, C for case 2.1. As shown in Figure 27, only a small 
portion of the image is detected as clutter, hence demonstrating that case 2.1 is 
substantially clutter-free. 

20 Figure 28 shows the ED frame for case 2.2. Note the clutter region in the 

vicinity of the transducer (mainly in rows 40 to 100). Figure 29 shows an image of M 
for case 2.2. As shown in Figure 29, the averaging smears most of the image, but the 
clutter region, residing mostly at the top of the image, remains relatively sharp. Figure 
30 shows an image of O for case 2.2. Note that pixels with high intensity indicate 

25 substantial changes in the gray-level during the cardiac cycle. Figures 31-32 show the 
gray-level histograms H M (Figure 31) and H 0 (Figure 32) for case 2.2. Five curves are 
shown in Figures 31-32, describing the original histogram (solid), the fitted curve 
(dashed), and the three calculated Gaussians (dotted). The intersection points between 
the Gaussians are 38 and 59 for H M , and 33 and 66 for H 0 . Thus, T M = 59 and 

30 7b = 33. As shown in Figures 31-32 the fitted curves are similar to the original 
histograms, and the three Gaussians are well separated. Figure 33 shows an image 
obtained by a pixel by pixel multiplication of M by the clutter binary matrix, C for 
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case 2.2. As shown in Figures 33, most of the clutter pixels are recognized, whereas 
most of the non-clutter pixels are not detected. 

Figure 34 shows the ED frame for case 2.13. Note the clutter region in the 
vicinity of the transducer (rows 1 to 60) and at the bottom left (mainly in rows 400 to 

5 450, columns 150 to 280). Figure 35 shows an image of M for case 2.13. As shown in 
Figure 35, the averaging smears most of the image, but the clutter region, residing 
mostly at the top and bottom left, remains relatively sharp. Figure 36 shows an image 
obtained by a pixel by pixel multiplication of M by the clutter binary matrix, C for 
case 2.13. As shown in Figures 36, most of the clutter pixels are recognized, whereas 

1 0 most of the non-clutter pixels are not detected. 



Discussion 

The present embodiments automatically detect clutter pixels, with high level of 
efficiency, as quantitatively demonstrated using the values of the ratios P d and P fat 

15 (mean P d value of 91.3 ± 13.6 %, and mean P fat value of 0.14 ± 0.24 %). Visual 
examination of all 16 test cases, and meticulous examination of the three specific 
reference cases further support the ability of the present embodiments to automatically 
detect clutter. In cases 2.2 and 2.13, which are affected by clutter, the majority of the 
clutter pixels were correctly detected, with minimal (or no) false detection. 

20 Inspection of the matrices M and O brings additional insight into the 

procedure. Clutter regions (detected visually) are characterized by low O values and 
high M values (otherwise, the effect of the clutter is negligible). Moreover, unlike 
other regions of the image, the clutter regions in M still appear grainy, despite the 
time-averaging. This is explained by the fact that clutter is characterized by small 

25 variations. This can be used as an additional criterion for clutter detection. 

Furthermore, assessment of both H M and Ho shows that the adaptive method 
for selecting threshold estimation provides a good fit to the histogram with well- 
separated Gaussians. 



30 EXAMPLES 

Left Ventricular Outline Detection 

The teachings of the present embodiments were employed for outlining the left 
ventricle. As further demonstrated below, the application of the teachings of the 
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present embodiments allows an automatic outlining of the left ventricle to a high level 
of accuracy. 

Materials and Methods 

5 Image Processing 

According to various exemplary embodiments of the present invention, the 
detection of the left ventricular outlines can be done by executing several method 
steps. Reference is now made to Figure 37, which is a flowchart diagram of a method 
of outlining a region, according to a preferred embodiment of the present invention. A 

10 description of the specific steps employed in the present example accompanies the 
general description of the method. As stated hereinabove, the method steps described 
hereinbelow can be executed either contemporaneously or sequentially in many 
combinations or orders of execution, and the ordering of the flowchart of Figure 37 is 
not intended to be limiting. 

15 Hence, the method preferably begins at step 100 in which a set-average 

intensity value is calculated for each picture-element (a pixel in the present example) 
to provide the average intensity matrix, M. The method preferably continues to step 
101 in which Hm, the intensity histogram of M, is fitted to a sum of a plurality of 
localized functions, which are used to define the thresholds T\ and T2 as intersection 

20 points between the localized functions. 

In the present example, the three Gaussians of Examples 1 and 2 (see Equation 
9) were used. The lower threshold, T\ 9 was set to the intersection point between G\ 
and G 2 , and the higher threshold, T 2 , was set to the intersection point between Gi and 
G 3 . 

25 Optionally and preferably, the method continues to steps 102 and/or 103 in 

which picture-elements corresponding to clutter (step 102) or the papillary muscles 
(step 103) are removed or interpolated. The removal of clutter can be done in any way 
known in the art. Alternatively, the clutter can be removed using apparatus 10 or the 
clutter removal method as exemplified in Figure 1 9. 

30 The papillary muscles are preferably detected under the postulation that they 

appear as small closed high-intensity regions within the left ventricular cavity. The 
assumption of high-intensity leads to the capability of locating the papillary muscles in 
M, rather than looking for them in each frame separately. 
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Following is a preferred procedure for detecting the papillary muscles. Hence, 
according to the presently preferred embodiment of the invention the thresholds T\ and 
T 2 are applied to M, in a Boolean manner to obtain two binary matrices B\ and B 2 \ 

BM=[M{m,n)>T,] 
B 2 (m,n)=[M( m ,n)>T 2 ]. 

Figures 38a-b are representative examples of images of the binary matrices B\ (Figure 
38a) and B 2 (Figure 38b). Once the binary matrices are constructed a series of 
morphological filters, such as, but not limited to, closing and opening are applied so as 
to smooth and fill small gaps in B\ and B 2 . Referring to the examples in Figure 38a-b, 
the application of the morphological filters results in a disappearance of the small 
closed regions in Figure 38a, and all but two closed regions in Figure 38b. 

The binary matrices B\ and B 2 are preferably used for calculating two 
additional binary matrices, Pi and P 2 , in which 1 's preferably define picture-elements 
belonging to the papillary muscles, and O's preferably define other picture-elements. 
The calculation of Pi and P 2 is preferably as follows: 

For/ = 1,2: 

(1) Randomly selecting an element (m p ,n p ) in B it whose value is 1 . 

(2) Defining a matrix R of the same size as B\, in which the values for all 
the elements but (m p ,n p ) are O's, while the value for (m p ,n p ) is 1. 

(3) Locating all the l's in B\ that belong to the closed region (of l's) in 
which (m p ,n p ) resides, and recording this region in R. This step can be done, for 
example, using a region-growing procedure. 

(4) Subtracting R from B„ thereby ensuring that the same region is not 
selected more than once, and, if the sum of R{m,n) is lower than a threshold R max , 
setting P\{m,n) = P\{m,n) + R(m,n). In the present example, a value of 500 was 
selected for R max . 

(6) Repeating steps (1) to (4), until Bi(m,n) contains only O's. 

Pi and P 2 are preferably merged to obtain a binary matrix P, defining the 
location of the papillary muscles: 

P(m,n) = P l (m i n)vP 2 (m,n) (EQ. 20) 

Figure 39 is a representative example of an image of the binary matrix P. As 
shown in Figure 39, only two regions of non-zero matrix-elements are left in P. The 
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bottom left region resides outside the left ventricle and the upper middle region is 
identified as corresponding to the papillary muscles. 

Referring again to the flowchart of Figure 37, in step 104 an origin (m/,«,) of 
the cine loop is defined and used for transforming the grid onto which the images are 
defined into polar representation. In the present example, the origin was defined as the 
center of mass of the left ventricle. 

Following is a preferred procedure for defining (m h n,): 

(1) Transform each element (m,n) in the imaging plane into a polar 
representation (/%#), using a preliminary reference point (m ri n r ). This results in two 
matrices, r and G, defined as follows: 



r(m,n)= ^j(m-m r f +(n-n r ) 2 

0(m, n) = arctarJ — — — 
ym-m r ) 



(EQ.21) 



The values in the r and 6 matrices can be further digitized to the required axial and 
lateral resolution, which can be the size of each range-gate and angular-gate. The 
digitized matrices are denoted and 

(2) Rejecting clutter regions and the papillary muscles, thereby obtaining a 
"free" matrix M/m,n): 

M / (m i n)=M(m,n) (\-P(m i n))i\-C(m i n)) (EQ. 22) 

where P is defined above (see Equation 20) and C is defined in Example 2 (see 
Equation 17). 

(3) Applying T\ to M/ f thereby to obtain a binary matrix B u and applying 
one or more morphological filters to B\, for smoothing and filling gaps within B\. 

(4) Constructing a boundary matrix in which the boundaries in B\, and, 
optionally and preferably the boundaries of the imaged region are stored. This can be 
done, for example, using the morphological filter "Erode", as follows: 



0(m, n) = B l (m, n) - Erode 



B x (m,n\ 



(EQ. 23) 



with a similar expression for the boundaries of the imaged region. 



49 

(5) Calculating an array /(ft/), defining the distance from the preliminary 
reference point to the inner boundary. Specifically, for each angular-gate 0 d , /(ft/) is 
calculated as follows: 



where min nz represents a function which returns the minimal non-zero value of a 
matrix. 

(6) Calculating the Cartesian location (m c ,n c ) of each point along the inner 
boundary defined by /(ft/). 

(7) Defining m, as the average value of m c and «, as the average value of n c . 
Once the origin («*,«/) is defined, the matrices r and 0 can be redefined using 

Equation 21 in which the preliminary reference point {m r ,n r ) is replaced by the 
calculated origin (/»,,«/). 

According to a preferred embodiment of the present invention the method 
continues to step 105 in which a region of the image in which the cusps of the mitral 
valve reside. In other words, in step 105 a region-of-interest containing the mitral 
valve is detected prior to the actual detection thereof. The advantage of this step is 
that it reduces the probability of falsely identifying the mitral valve. 

The detection of the mitral valve region-of-interest is based on the hypothesis 
that the location of the cusps of mitral valve, as seen within the imaging plane, varies 
dramatically during the cardiac cycle. As a result, the relevant pixels (in each frame) 
are expected to show large deviations from the mean intensity over time, M. In 
addition, these pixels also tend to have high pixel intensities. On the other hand, the 
above definition also applies to the edges of the cardiac muscle. According to a 
preferred embodiment of the present invention this fact is addressed in two phases. In 
a first phase maximal dimensions for the cusps are assumed, and in a second phase 
additional filters are applied for the final determination. 

Hence, the initial region-of-interest determination is performed by: 

(1) Defining an initial matrix Vj(m,ri) = 0V/w,«. 
For each p: 

(2) Calculating D(m,n), the relative deviation of each pixel from M(m,n): 




(EQ. 24) 



D(m,n) = 



\F p (m,n)-M(m,n 



(EQ. 25) 
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(3) Setting elements satisfying D(m,n) < T 2 to 0. 

(4) Applying a constant threshold, T D , to D(m,n), to obtain £>t(/w,h): 

D T (m, n) = (Z)(m, n)>T D ) (EQ. 26) 

in the present example a value of 1.0 was selected for T D . 

(5) Constructing a binary matrix //, defining high-intensity regions in 
F p (m,n), based on T\ and r 2 : 

H(m,n)={F p (m,n)>\(T x +T 2 )) (EQ. 27) 

(6) Updating the matrix V(m,ri), by adding 1 to pixels (m p ,n p ), for which all 
the following conditions apply: 

(a) D 7 (m p ,n p ) = \; and 

(b) when examining a block HBim^rn,) of H(m,n), of a constant size 
(H\ x // 2 ), centered at (m p ,n p ), one of the following statements is true: 



£ Y,H B {m b ,n h ) 



1 "2 



n b =H 2 -2 



-o 



(EQ. 28) 



(EQ. 29) 



each bracketed expression in Equations 28-29 refers to a different corner of H B . The 
above expressions mean that at least two opposite corners of the block show relatively 



low intensity, such that the local width does not exceed tJh* + H 2 2 . This defines the 
maximal allowable width of the mitral valve cusps (used for discerning the cusps from 
the cardiac muscle). 

Once all frames are processed, V t is transformed into a binary matrix: 

^(m 9 n)=(K i (m,«)>l), (EQ. 30) 

and its largest closed region is defined as the mitral valve region-of-interest, V. 

This can be done, for example, by employing a procedure which is similar to 
the procedure employed for obtaining the matrices P\ and P 2 . Hence, according to a 
preferred embodiment of the present invention, V is defined by: (i) randomly selecting 
an element (m p ,n p ) in V h whose value is 1 ; (ii) defining a matrix R, in which the values 
for all the elements but (m p ,n p ) are 0's, and the value for (m p ,n p ) is 1; using a region 
growing procedure for storing in R all the 1 's in F, which belong to the closed region 
of l's in which (m p ,n p ) resides; subtract R from K„ so as not to select the same region 
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more than once; and (iii) if the sum of R exceeds a threshold V s , then setting V— R and 
updating V s to be the sum of /?, where the initial value of V s can be set to zero. The 
process is iteratively repeated until V\{m,ri) = 0 for all (m,ri). 

The above process is demonstrated in Figures 40a-b, showing an image of M 
5 (Figure 40a) and the mitral valve region-of-interest (Figure 40b). Comparing Figure 
40a and Figure 40b, one can recognize that the region-of-interest shown in Figure 40b 
encloses the mitral valve. 

The region-of-interest, once determined, is preferably transformed into polar 
representation. For apical echocardiograph images, for example, the mitral valve is 
10 typically located at the bottom of the left ventricle. Thus, the transformation is 
performed by (i) locating in V the first and the last columns having non-zero elements; 
(ii) locating the highest row-index of non-zero elements in these columns; and 
calculating the following angular-gates using the matrix 6d(m,n): 

f d ) f f J (EQ. 31) 

1 5 where N/ and Ni are the aforementioned first and last columns, and Mr and M/ are their 
highest row-index of non-zero elements, respectively. The angular gates ^and Oi thus 
span the region-of-interest of the mitral valve. 

Whether or not step 105 is employed, the method preferably continues to step 
106 in which the inner and/or outer boundaries of the left ventricle are detected. The 
20 boundaries are preferably determined in the polar representation of the grid and 
subsequently, if desired, are transformed to the original representation. The inner and 
outer boundaries of the pth image are denoted hereinbelow byR' n (O d ,p) and 
R out (6 d ,p) respectively. 

The determination of boundaries is preferably preceded by the removal of the 
25 clutter and the papillary muscles are from each frame, thereby to provide a filtered 
matrix Ff'(m,n). This is preferably done by setting all the clutter pixels and the 
papillary muscles' regions to 0: 

Ff (m, n) = F p {m, b). (l - P(m 9 «)) ■ (l - C(m s b)) . (EQ. 32) 

Optionally and preferably, determination of boundaries is also preceded by a 
30 procedure in which the mitral valve cusps are located and removed from each filtered 
matrix. This can be done similarly to the procure of detecting the mitral valve region- 
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of-interest 5 with the exceptions that T\ is used as a high-pass filter for the relative 
deviation matrix, D, and V is used as the region-of-interest in which the cusps are 
located. 

Thus, according to the presently preferred embodiment of the invention, the 
5 cusps are located by: (i) calculating D according to Equation 25 while setting elements 
satisfying D(m,ri) <T\ to 0; (ii) constructing a matrix V p defined as 
V p = (D(m,n) > T D ), where T D is the constant defined in Equation 26 above; (iii) 
multiplying V p {m,ri) by V{m,n) which, as stated, stores the mitral valve region-of- 
interest, thereby substantially preventing false identification; (iv) multiplying 
10 Ff'(m s n)by \-V p (m 9 n), thereby removing the mitral valve from the filtered matrix; 

(v) locating in V p the first and the last columns having non-zero elements; (vi) locating 
in these columns the highest row-index of non-zero elements; and (vii) selecting the 
(rectangular) region defined by the located non-zero elements, thereby locating the 
region of the mitral valve cusps. 
15 The radial-gate and angular-gate of the mitral valve cusps can be calculated 

using the radial and angular matrices. For example, if the digitized matrices, r</ and @d, 
are used, radial- and angular-gates of the mitral valve cusps are preferably calculated 
as follows: 

0?=O d {M P f>X P f) (eq.33) 

20 where ( M p f , N£ ) and ( Mf , Nf ) define the corners of the mitral valve cusps region of 
the /?th frame. 

The process of determining the boundaries preferably begins in a construction 
of one or two binary matrices, Bf (m, n) and B£(m,n), for each frame, p. This is 
preferably done by an application of a strong low-pass filter, followed by a Boolean 
25 application of the thresholds T\ and T 2 to Ff"(m,n). The process continues by 

locating the inner and/or outer boundaries of the binary matrices Bf(m,n) and 
B%{m,n) of the pth. frame. Generally, this is done by traversing rays emanating from 
the origin (/»,-,«,), at every angle. 
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Figure 41a is a schematic illustration of a preferred process for determining the 
left ventricle boundaries, based on the binary matrices ovBf{m,n) and B[(m y n). The 
left ventricular cardiac muscle is shown in Figure 41a as a gray region and designated 
by numeral 120. Rays 122 emanating from the origin {m h n t ) are marked as dashed 
5 lines. The intersection points of the rays with the inner and outer boundary are 
designated by numerals 126 and 124, respectively. Note that since the mitral valve is 
removed from the image, no intersection points are found within the data region for 
the corresponding rays (e.g., ray 122'). 

According to a preferred embodiment of the present invention, the locations of 

10 both the inner boundary and the outer boundary are described by their local distance 
from the origin, along each ray. For each ray, the range to the inner boundary is 
preferably defined as the range to the first non-zero value along the ray, and the range 
to the outer boundary is preferably defined as the range to the first zero-pixel beyond 
the inner boundary. As will be appreciated by the skilled artisan, in some regions, no 

15 boundary is found within the data region. This may be the case, for example, in the 
region of the mitral valve, which, as stated, is preferably removed prior to searching 
for the boundaries. Such regions are preferably marked, e.g., by assigning the 
corresponding matrix element with a value other than 0 or 1 (in the present example a 
value of -1 was selected). 

20 Following is a detailed description of a preferred process for locating the inner 

and outer boundaries of the binary matrices, which is preferably performed for each 
frame p and for each / = 1, 2. 

(1) Smoothing and filling gaps m.B?{pt,n) by applying a set of 
Morphological filters thereto. 

25 (2) Constructing a boundary matrix in which the boundaries in B h and, 

optionally and preferably the boundaries of the imaged region are stored. This can be 
done, for example, using the morphological filter "Erode", as follows: 

1 1 1 
1 1 1 



0(m, n) = B i (m, n) - Erode 



1 1 1 
1 1 1 
1 1 1 



(EQ. 34) 



with a similar expression for the boundaries of the imaged region. 
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(3) Calculating a matrix Rf (O d , p) , defining the distance from the origin to 
the inner boundary. Specifically, for each angular-gate 0 d , R i "($ di p) is calculated as 
follows: 

^(^,/?)=min W2 {r rf (m,«).(^(m,«)=^) 0(m, W )} (EQ. 35) 

where, min nz returns a minimal non-zero element, or -1 if no non-zero element is 
found. 

(4) Calculating a matrix R°"(0 a ,p), defining the distance from the origin 
to the outer boundary. Specifically, for each 0 A R° ut (0 d ,p) is calculated as follows: 

^rfe,p) = nmin„{r rf (i W> 4(^(m,/i) = 6>J.O(jfi,w)} (EQ. 36) 

where nmin nz returns the next-to-minimal non-zero element, or -1 if no significant 
result can be obtained (e.g., when there is no such element, or the next-to-minimal 
value is substantially close to the minimal value). 

Once the matrices R^{9 d9 p)^(0 di p) 9 R^(0 d9 p) and R? ! (0 di p) are formed, 
the process preferably combines them so as to obtain the boundaries of the left 
ventricle. While reducing the present invention to practice, it was found by the present 
Inventors that inner boundary, referred to herein as R in (0 d ,p), typically coincides with 
the boundary delineated by R[ n (e di p). Further it was found by the present Inventors 
that if the difference between R™'(0 d ,p) and R?'(0 d ,p) is sufficiently large, the 
outer boundary, referred to herein as R oul (0 d , p) , equals the smaller value 
between Rf"(0 a ,p) andR?'(0 d ,p). On the other hand, for angular-gates in which 
R r'(&d>p) and R^(0 d ,p) are close, R out (0 d> p) equals the larger value therebetween. 

Thus, according to a preferred embodiment of the present invention, 
R in {0 <J ,p)=Ri n (0 di p),and 



"™ \min{Rr(0 d ,p\Rr{O d ,p)} otherwh 



< R 

- mozdif (eq. 37) 
otherwise 



where, R maxdif is the maximal difference between R° ul {e d ,p) and Rl ut {6 di p) i 
™**%fa>p) = \RT*(0 a ,p)-R?{0 d9 p\ In the present example, the-1 values in 
both R? ul (0 d ,p) and R^'fa^p) were replaced by the maximal possible values. 
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Additionally, the elements in JT'fo,/?), which equal the maximal possible value are 
preferably set to their original value, -1 . 

Once R' n and R out are calculated, the method preferably continues to step 107 in 
which a set of temporal and spatial filters are applied so as to reject non-valid local 
information in R in (0 d ,p) and R out (0 d ,p). The rejected information is preferably 
replaced with interpolated data. 

With respect to the inner boundary, the temporal and spatial filters are 
preferably applied according to the following procedures: 

(1) Detecting angles Q d , for which the standard deviation in R in (0 d ,p) 
exceeds a threshold cr max , and replacing all the R ,n (0 d ,p) values for the corresponding 
angles by -1. cj max is typically set to the median value of the standard deviations, 
multiplied by a constant, o 0 . In the present example, a value of 3.5 was selected for 
a 0 . 

(2) Replacing all the elements belonging to the mitral valve (elements 
characterized by 0 p f < 0 d < Qj> ) by -1 . 

(3) Applying a strong two-dimensional low-pass filter toR in (O d> p), 
thereby obtaining^ (9 d ,p), and replacing elements, for which the ratio | R'"(O d ,p) - 
r u>f(0<i>p) I / R'Hpf(^^p) is larger than a threshold A max , by -1. A typical value for 
A max is from 0.05 to 0.2. In the present example, a value of 0.15 was selected for A max . 

(4) Performing two-dimensional interpolation and replacing elements in 
R m (0 d ,p) having a value of -1, with interpolated values. Preferably, a cyclic 
interpolation is performed so as to coincide with the cyclic nature of the data. 
Specifically, if an integral number of cardiac cycles have been used, the interpolation 
is such that R m (0 d ,p) is cyclic in both axes, otherwise, R in (9 d ,p) is cyclic in 9 d . This 
step is preferably supplemented by a two-dimensional smoothing procedure (e.g., via 
spline fitting) so as to smooth R in (0 d ,p). 

(5) Connecting the points (r/,0>) and (r t ' 9 0f), which, as stated, define in 
polar representation the base of the mitral valve cusps, by a straight line. 
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With respect to the outer boundary, steps (l)-(5) of the inner boundary are 
preferably repeated, for a matrix representing the local muscle width, defined as 
W(O d ,p) = R ou \O d ,p) - R [n (0 d ,p). Once steps (l)-(5) are repeated, R out (0 d ,p) is 
preferably reconstructed according to the formula: 

R°»(e d ,p) = R'"{9 d ,p)+W(e d ,p). (EQ. 38) 

According to a preferred embodiment of the present invention the outer 
boundary is limited by the extensions of the straight line connecting the base of the 
mitral valve cusps. This can be done by replacing outer boundaries exceeding the 
straight line with a continuation of the line. 

Figure 41b is a schematic illustration of left ventricle after the connection of 
the base of the mitral valve cusps 130 by a straight line 134. Also shown in Figure 
41b are extensions 128 of line into the cardiac muscle, and the left atrium 132. As will 
be appreciated by one ordinarily skilled in the art, the process of outlining left atrium 
132 is similar to the present process in which left atrium 132 was outlined. 

In optional and preferred step 108, the boundaries are transformed back from 
polar representation to the original or any other representation. For example, for the 
case of a Cartesian coordinate system, the bounties can be transformed by: 

( 1 ) Calculating a binary matrix I(m,ri) defined as: 

/(m,«)==X(o(^«)^^fe^)>fe(^«) = ^)- (EQ. 39) 



20 



25 



The matrix / thus includes l's for elements within the inner boundary. Note that in 
Equation 39, the symbol 0 d is used both as an index 6 d and as a matrix 6<j(m,ri). 

(2) Determining the inner boundary in Cartesian coordinates, R'"(m, n,p), 
using the morphological filter "Erode": 



R'"(m,n,p)= l(m,n)- Erode 



(EQ. 40) 



(3) Repeating steps (1) and (2) for the R ou \ thereby obtaining the outer 
boundary in Cartesian coordinate R^fa, n,p). 
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Data acquisition 

Seven cine-loops of echocardiograph images recorded of three different 
patients were selected. The cine-loops were in apical four-chamber and apical two- 
chamber views. The data has been collected using a Vivid™ 3 imaging platform, by 
5 GE Medical Systems Ultrasound. 

Results 

Figures 42a-e and 43a-e show results for two representative cine-loops 
collected at a high frame rate of 75 frames per second: case 3.1 (Figures 42a-e) and 

10 case 3.2 (Figures 43a-e). Case 3.1 contained a single cardiac cycle at a heart rate of 65 
beats/min, and case 3.2 contained a single cardiac cycle at a heart rate of 71 beats/min. 
Case 3.1 was formed of 50 frames and case 3.2 was formed of 42 frames. 

Shown in Figures 42a-43e are epicardial boundaries 142 and endocardial 
boundaries 144 and frame Nos. 1 (Figures 42a and 43a), 1 1 (Figures 42b and 43b), 21 

15 (Figures 42c and 43c), 31 (Figures 42d and 43d), and 41 (Figures 42e and 43e). In 
each frame, endocardial boundary 144 is marked as an inner red closed line, and 
epicardial boundary 142 is marked as an outer magenta closed line. As shown, the 
results are in complete agreement with the visually determined outlines. In regions 
where the cardiac muscle extends outside of the scanned data region, the outlines 

20 follow the boundaries of the data region 

Discussion 

The present embodiments successfully detect the outlines of the cardiac muscle 
within the imaging plane of echocardiograph images. Both temporal and spatial filters 

25 have been used so as to handle low signal-to-noise ratio, which is a common drawback 
in echocardiography. Special treatment was provided for clutter, so as to allow 
outlining to go through clutter regions. The cusps of the mitral valve and the papillary 
muscles have been automatically detected and removed from the images prior to the 
outlining process. The present embodiments have a minimal number of assumptions 

30 regarding the nature of the cardiac muscle, and hence support both intra-patient and 
inter-patient variations. 

Careful examination of the results reveals that the detected outlines closely 
match the visual estimation for the location of the cardiac muscle in all the cine-loops 
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in the test-group. The present embodiments successfully handle clutter, as well as the 
complex motion of the mitral valve and the papillary muscles. Small variations appear 
mainly in the vicinity of the mitral valve. This occurs when a small part of the cusps 
is not rejected by the pre-filtering, and is thus marked as a part of the cardiac muscle 
itself. 

The present embodiments provide an automated estimate for the location of 
both the endocardial and the epicardial boundaries of the left ventricle. The quality of 
the results is sufficient for the extraction of various local quantitative functional 
parameters, which are indicators for muscle contractility and viability. 

It is appreciated that certain features of the invention, which are, for clarity, 
described in the context of separate embodiments, may also be provided in 
combination in a single embodiment. Conversely, various features of the invention, 
which are, for brevity, described in the context of a single embodiment, may also be 
provided separately or in any suitable subcombination. 

Although the invention has been described in conjunction with specific 
embodiments thereof, it is evident that many alternatives, modifications and variations 
will be apparent to those skilled in the art. Accordingly, it is intended to embrace all 
such alternatives, modifications and variations that fall within the spirit and broad 
scope of the appended claims. All publications, patents and patent applications 
mentioned in this specification are herein incorporated in their entirety by reference 
into the specification, to the same extent as if each individual publication, patent or 
patent application was specifically and individually indicated to be incorporated herein 
by reference. In addition, citation or identification of any reference in this application 
shall not be construed as an admission that such reference is available as prior art to 
the present invention. 



